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SECTION  I 


INTRODUCTION 


The  purpose  of  this  study  was  to  determine  the  real-time  hardware 
requirements  for  computer  generation  of  nonlinear  textured  terrain  imagery. 
The  approach  is  based  on  algorithms  and  software  developed  by  Honeywell 
under  an  Independent  Research  and  Development  (IR&D)  program.  This 
software  was  implemented  on  the  Air  Force  Human  Resources  Laboratory 
(AFHRL)  Simulation  Training  and  Advanced  Research  System  (STARS) 
computer  facility  for  nonreal-time  generation  of  imagery  during  Phase  I 
of  this  program.  The  quality  of  the  resulting  imagery  was  considered 
sufficient  to  warrant  an  investigation  of  real-time  hardware  requirements 
during  Phase  II.  This  report  discusses  that  investigation. 

This  effort  is  part  of  a  larger  AFHRL  program  to  apply  computer  image 
generation  (CIG)  to  techniques  for  use  in  visual  and  electro-optical  (E/O) 
sensor  simulation  systems  for  trainer  systems.  Emphasis  on  low-level 
flying  and  the  need  to  incorporate  training  for  target  detection  and  identi¬ 
fication  tasks  places  new  requirements  for  realism  in  these  simulations. 
Previous  studies  have  established  the  applicability  of  CIG  technology  to 
visual  and  E/O  sensor  simulation  and  have  developed  algorithms  for  the 
simulation  of  atmospheric  effects  and  sensor  transfer  characteristics. 

These  studies  were  based  on  existing  CIG  techniques  which  use  polygonal 
representation  of  surfaces.  Smooth  shading  techniques  were  incorporated 
for  added  realism. 

Polygon  methods  are  not  readily  extended  to  include  texturing  of  surfaces. 
Furthermore,  the  computational  requirements  for  polygonal  terrain 
representation  would  be  considerable.  Although  smooth  shading  can  reduce 
or  eliminate  facetting  effects  in  the  interior  of  objects,  silhouettes  and 
boundaries  are  still  composed  of  straight  line  segments  which  detract  from 
the  realism  of  the  simulation. 

The  alternative  approach  which  is  the  basis  for  this  study  represents 
surfaces  as  collections  of  curvilinear,  bivariate  patches.  In  general,  any 
surface  in  three-dimensional  space  can  be  represented  as  a  function  of  two 
variables.  For  example,  a  terrain  surface  is  determined  by  elevation 
above  some  datum.  The  two  independent  variables  represent  coordinates 
in  the  datum,  or  ground  plane.  Other  bivariate  functions,  representing 
surface  material,  texture,  etc.  ,  can  be  defined  with  respect  to  the  same 
coordinates,  thus  allowing  texturing  of  surfaces. 
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The  particular  representation  used  here  for  terrain  is  the  bicubic  surface 
patch,  a  polynomial  approach  which  has  recently  received  extensive  develop¬ 
ment  for  use  in  computer-aided  design  of  surfaces.  The  bicubic  patch 
representation  has  first  or  second  derivative  continuity  across  patch 
boundaries,  depending  on  the  particular  bicubic  representation  used  (see 
Appendix  A).  Consequently,  smooth  shading  and  curved  silhouettes  and 
boundaries  are  inherent  in  the  approach. 

The  generation  of  a  perspective  view  of  a  surface  consisting  of  bicubic 
patches,  as  required  for  visual  and  E/O  sensor  simulation,  requires 
different  methods  than  those  used  for  polygon  surfaces.  The  method  used 
here  is  the  bicubic  subdivision  method  of  E.  Catmull.  ^  Catmull's  method 
is  a  fast,  recursive  subdivision  of  the  patch  into  smaller  and  smaller  sub¬ 
patches,  until  each  subpatch  is  small  enough  to  be  displayed.  Because  the 
order  in  which  the  subpatches  are  ready  for  display  is  random,  the  usual 
scanline  approach  in  which  the  imagery  is  produced  a  line  at  a  time  in 
synchronism  with  the  raster  scan  of  the  display  is  not  applicable.  Rather, 
the  image  is  assembled  in  a  special  memory  called  a  frame  buffer. 

The  frame  buffer  contains  the  entire  image  at  display  time  and  is  scanned 
out  in  raster  sequence.  Catmull  extended  the  frame  buffer  concept,  originally 
proposed  by  M.  Newell,  l2-*  to  include  the  distance  from  the  viewpoint  to  each 
displayed  point.  The  extended  frame  buffer  is  called  a  Z-buffer,  which 
provides  a  simple  solution  to  the  hidden  surface  problem  by  resolving 
conflicts  based  on  distance  from  the  viewpoint. 

The  Z-buffer  also  includes  the  following  advantages.  Computing  requirements 
increase  only  linearly  with  scene  complexity.  For  other  priority  algorithms, 
computing  costs  increase  faster  than  linearly  with,  say,  the  number  of  edges 
required  for  the  scene.  Also,  different  and  more  suitable  object  represen¬ 
tations  may  be  used  within  the  same  system.  For  example,  bicubic  patches 
may  be  used  for  terrain,  polygons  for  buildings,  and  quadratics  for  storage 
tanks  and  water  towers.  Since  distance  to  each  picture  element  is  available, 
atmospheric  propagation  effects  can  be  included  by  table  look-up  at  display 


surfaces.  UTEC-CSc-74-133,  University  of  Utah,  December  1974. 


r2l 

JNewell,  M.  E.  ,  Newell,  R.G.  and  Sancha,  T.  L.  "A  Solution  to  the 
Hidden  Surface  Problem,  "  Proc.  ACM  Nat.  Conf.  ,  August  1  972. 


time.  Finally,  since  the  neighbors  of  each  picture  element  are  available, 
two-dimensional  interpolation  and  a  form  of  anti-aliasing  filtering  are 
easily  incorporated. 

In  order  to  study  the  Honeywell  approach  to  CIG,  AFHRL  at  Wright- 
Patterson  AFB  contracted  with  the  Honeywell  Systems  and  Research  Center 
for  a  software  implementation  of  the  algorithm  for  evaluation.  Under  Phase  I 
of  contract  No.  F33615-77-C-1173,  "A  Study  of  Real-Time  Feasibility  for 
Generation  of  Nonlinear  Textured  Terrain,  "  Honeywell  transferred  a  working 
version  of  its  CIG  software  and  data  to  the  AFHRL  facility  at  Wright-Patterson 
AFB,  Ohio.  Under  Phase  II  of  this  contract,  Honeywell  defined  the  real-time 
hardware  organization  required  to  implement  the  algorithm.  The  study 
established  the  requirements  using  a  combination  of  analytical  methods  and 
computer  simulation,  analyzed  several  system  configurations,  and  performed 
a  detailed  study  of  memory  and  computational  devices  required.  This  report 
covers  the  technical  effort  for  Phase  II  of  the  contract. 
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SECTION  II 


PROBLEM  AND  APPROACH 


The  Honeywell  CIG  approach  combines  several  innovative  ideas.  These 
include  using  a  frame  buffer  and  depth  buffer  to  assemble  the  components 
of  a  scene,  using  curved  surfaces  instead  of  polygons  for  terrain  represen¬ 
tation,  applying  texture  patterns  to  the  surfaces,  and  using  polygons  to 
represent  cultural  objects  in  the  same  scene. 

The  basic  elements  of  the  CIG  problem  for  terrain  are  shown  in  Figure  1, 
They  consist  of  a  viewpoint  represented  by  an  eye,  a  "window,  "  the  terrain, 
and  a  line  of  sight  (LOS).  The  window  may  be  an  actual  window  in  the  case 
of  visual  simulation,  or  the  focal  plane  in  the  case  of  a  camera  such  as 
forward-looking  infrared  (FLIR)  or  low-light-level  television  (LLLTV) 
camera. 

In  simulator  visual  systems,  the  window  is  replaced  by  a  cathode  ray  tube 
(CRT)  or  similar  display  device  with  suitable  optics  to  provide  a  virtual 
infinity  image  at  the  viewpoint.  The  scene  viewed  through  the  window  is 
synthesized  by  computer  from  a  data  base.  The  data  base  consists  of 
terrain  elevation  samples  and  descriptors  of  the  material  at  the  surface 
of  the  ground  as  a  function  of  geographical  coordinates. 


Figure  1.  Terrain  Image  Generation  Problem 


The  synthesized  image  consists  of  a  rectangular  array,  or  raster,  of  intensity 
samples.  The  intensity  or  color  displayed  at  a  particular  raster  element  is 
determined  by  the  material  characteristics  (reflectance,  emittanee,  temper¬ 
ature,  etc.)  of  the  first  surface  intersected  by  the  LOS  from  the  viewpoint 
through  the  raster  element.  It  is  also  determined  by  illumination,  surface 
orientation,  and  attenuation  along  the  LOS.  For  FL1R  and  LLLTY,  the 
displayed  intensity  depends  also  on  camera  resolution,  sensitivity,  and 
noise  characteristics. 

The  raster  element  should  not  be  considered  as  a  point,  however,  but  as 
a  small  element  of  area  on  the  display.  The  solid  angle  determined  by  the 
viewpoint  and  the  raster  element  will  intersect  the  terrain  over  some  patch 
of  terrain  (Figure  2).  The  patch  is  the  projection  of  the  raster  element 
onto  the  terrain  surface.  The  area  of  the  patch  will  vary  with  distance 
and  obliquity  of  the  surface. 

The  displayed  intensity  is  determined  by  the  average  brightness  over  the 
path  and  not  by  the  brightness  at  a  point,  as  suggested  by  Figure  1. 

Required,  then,  are  algorithms  for  sampling  for  average  brightness  of  the 
ground  patch  corresponding  to  each  raster  element  within  the  field  of  view. 

As  the  viewpoint  and  orientation  of  the  window  change,  the  projection  of 
ground  patches  on  the  display  window  will  change  accordingly.  The 
algorithms  must  be  computationally  simple  enough  to  allow  for  real-time 
tracking  of  these  changes  with  realistic  hardware  speeds. 
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BICUBIC  SPLINES 


Bicubic  splines  provide  a  smooth  polynomial  interpolation  of  the  terrain 
elevation  samples  in  the  Defense  Mapping  Agency  (DMA)  source  data. 

Each  spline  approximates  a  patch  of  terrain  overlaying  a  rectangular  area 
in  the  ground  plane  (Figure  3).  The  ground  plane  in  this  sense  is  an 
arbitrary  horizontal  datum  plane.  A  Cartesian  coordinate  system  is 
defined  with  the  Z  axis  aligned  vertically  and  the  X  and  Y  axes  aligned 
east  and  north,  respectively,  with  origin  directly  beneath  the  aircraft  in 
the  ground  plane.  The  terrain  patches  then  correspond  to  the  latitude/ 
longitude  grid  of  the  DMA  source  data. 

The  patch  boundaries  are  intersections  of  the  four  vertical  planes  along  the 
latitude /longitude  grid  lines  with  the  terrain  surface,  and  are  usually  curved 
due  to  the  bicubic  representation.  The  projections  of  these  curves  onto  the 
window  are  also  curved.  Consequently,  unlike  polygon  representations, 
the  four  vertices  cannot  simply  be  connected  with  straight  lines  in  the 
window  plane  to  determine  the  imaged  patch  boundaries.  The  approach 
taken  here  is  to  subdivide  the  patch  until  the  resultant  subpatches  each 
cover  one  picture  element. 
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A  bicubic  spline  provides  a  polynomial  representation  of  terrain  elevation 
z  of  third  degree  in  x  and  y  which  can  be  written  as 


2  2  2  3  3 

Z  =ao+alX  +  a2y  +V  +V  +Vy  +---+a15xy 


A  bicubic  has  16  arbitrary  coefficients  which  are  determined  by  a  minimum 
of  16  terrain  elevation  samples.  The  simplest  approach  uses  the  four 
corner  elevations  for  the  patch  and  the  neighboring  12  elevation  samples. 

The  bicubic  coefficients  can  be  simply  determined  by  arranging  the  16 
elevations  into  a  4  x  4  matrix  £  and  then  pre-multiplying  and  post-multiplying 
Z  by  a  4  x  4  matrix  of  constants,  M,  and  its  transpose,  respectively. 
Appendix  A  discusses  bicubic  spline  mathematical  properties  in  more  detail. 

Depending  on  the  choice  of  M,  a  bicubic  may  be  either  interpolating  or 
approximating.  An  interpolating  spline  surface,  which  has  first  derivative 
continuity,  is  one  which  fits  the  true  surface  at  the  four  corners  of  the  patch 
and  interpolates  the  surface  elsewhere.  An  approximating  spline,  which 
has  second  derivative  continuity,  approximates  the  surface  everywhere. 

In  either  case,  the  resulting  surface  is  smooth  with  at  least  first  derivative 
continuity  everywhere,  including  along  patch  boundaries.  Thus,  the  normal 
to  the  surface  is  continuous.  This  is  important  because  shading  of  surfaces 
depends  on  the  angles  between  the  surface  normal  and  lines  to  the  viewpoint 
and  illumination  source.  The  human  eye  is  very  sensitive  to  spatial  dis¬ 
continuities  in  shading,  even  in  the  first  derivative  of  shading;  therefore, 
discontinuities  in  the  surface  normal  may  be  discernible.  The  effect  might 
be  particularly  objectionable,  and  in  fact  may  provide  false  cues  to  the 
trainee,  in  real-time  visual  simulation,  if  the  discontinuities  appear  to 
jump  or  crawl  as  the  angles  change. 


Special  interpolation  techniques  have  been  developed  to  reduce  these  effects 
when  representing  curved  surfaces  by  polygons.  Examples  are  GouraudMl 
shading  and  PhongL  J  shading.  An  advantage  of  bicubic  spline  representation 
of  surved  surfaces  is  that  the  necessary  continuity  is  inherent  in  the  approach. 
No  smoothing  approximations  are  necessary. 


^Gouraud,  H.  "Computer  Display  of  Curved  Surfaces,  "  Dept.  Comp. 

Science,  U.  of  Utah,  Tech.  Rept.  UTEC-CSc-71 -1 1 3,  June  1971. 

^Phong,  Bui  Tuong,  "illumination  for  Computer-Generated  Images,  "  Dept, 
of  Comp.  Science,  U.  of  Utah,  Tech.  Rept.  UTEC-CSc-73-1 29,  July,  1973. 
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Bicubics  are  the  splines  of  lowest  degree  which  have  this  property.  Although 
biquadratic  splines  require  fewer  coefficients,  they  cannot  provide  the 
desired  degree  of  continuity.  This  is  the  principal  reason  bicubic  splines 
were  chosen  for  terrain  representation. 

GROUND  PLANE  INTENSITY  REPRESENTATION 

The  previous  subsection  dealt  with  the  mathematical  representation  of  terrain 
for  CIG.  The  terrain  can  best  be  represented  by  bicubic  spline  functions; 
these  are  bivariate  polynomials  determined  by  the  data  base  terrain  elevation 
samples.  The  bicubic  representation  was  chosen  for  its  smoothness  properties. 

The  other  half  of  the  DMA  data  base  consists  of  planimetric  features  which 
can  be  represented  as  a  map  in  the  ground  plane  (Figure  4).  Each  surface 
material  (concrete,  grass,  rock,  sand,  etc.  )  is  represented  by  a  code  as  a 
function  of  latitude /longitude  grid  coordinates.  From  this  data,  a  map  like 
that  of  Figure  4  can  be  constructed  to  the  resolution  level  of  the  data  base. 

Such  a  map  can  also  be  made  from  any  two-dimensional  intensity  or  color 
distribution  such  as  a  photograph  or  sensor  image.  Aerial  photographs 
and  random  patterns  were  used  in  the  Phase  I  software  implementation. 

It  is  possible  to  assign  an  intensity  (representing  reflectance,  emittance, 
or  other  material  characteristics)  to  each  surface  material.  Each  intensity 
can  be  mapped  from  the  ground  plane  to  the  image  plane  by  mapping  the 
curved  terrain  surface  onto  the  image  plane  (Figure  5).  The  mapped 
intensity  value  can  then  be  modified  for  illumination,  temperature, 
atmospheric  attenuation,  noise,  etc.  to  give  the  displayed  intensity. 

The  ground  plane  intensity  values  can  be  specified  to  any  desired  spatial 
resolution.  In  particular,  this  concept  allows  texture  to  be  specified  in 
the  ground  plane  for  each  terrain  material  and  mapped  onto  the  image 
plane  with  correct  perspective,  occlusion,  and  variable  resolution  with 
distance. 

Texture  in  images  can  be  described  as  a  quasi-periodic  or  random  variation 
of  intensities  which  gives  a  surface  a  particular  character.  The  surface  is 
then  perceived  as  tree-covered,  or  as  a  plowed  field,  etc.,  by  the  viewer. 
Perspectively  correct  texture  can  also  give  the  impression  of  distance, 
surface  slant,  and  motion.  Consequently,  texture  should  add  significantly 
to  realism  and  depth  and  motion  cues  in  CIG. 
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Figure  4.  Ground  Plane  Map  of  Planimetric  Features 
BICUBIC  SURFACE  MAPPING 

The  approach  described  here  for  CIG  using  bicubic  spline  terrain  represen¬ 
tation  is  based  on  Catmull's  bicubic  spline  subdivision  algorithm.  The 
algorithm  is  described  in  detail  in  Appendix  B. 

Figure  6  illustrates  a  mapping  of  a  patch  onto  the  image  plane  and  two 
successive  stages  of  subdivision.  The  original  patch  corners  are  denoted 
A,  B,  C,  D.  The  projection  of  the  patch  boundary  curves  are  shown 
connecting  these  points.  Catmull's  algorithm  provides  the  new  boundary 
curves  for  the  binary  subdivision  of  the  first  patch. 
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Figure  6.  Mapping  of  Ground  Plane  onto  Image  Plane 

Each  bicubic  patch  is  successively  subdivided  into  four  subpatches  until 
each  of  the  resulting  subpatches,  when  projected  onto  the  image  plane, 
occupies  the  area  of  one  raster  element.  The  visibility  of  each  such  sub¬ 
patch  is  determined  and  the  corresponding  intensity  displayed.  Because 
the  image  is  constructed  in  a  scattered  fashion,  a  frame  buffer  is  used  to 
assemble  the  picture  prior  to  display.  If  distance  to  and  intensity  of  each 
subpatch  is  stored  in  the  buffer  for  each  raster  element,  the  visibility 
determination  is  greatly  simplified.  This  concept  is  called  a  Z-buffer. 

Several  sequences  of  scenes  have  been  generated  to  demonstrate  the  capa¬ 
bilities  of  the  program.  Figure  7  shows  successive  views  of  one  corner 
of  a  data  base  including  some  buildings  on  the  terrain.  Figure  8  is  a 
series  of  views  of  a  more  mountainous  portion  of  the  data  base. 

The  data  base  environment  for  Figures  7  and  8  consists  of  a  hilly  area  of 
southwest  Nevada.  The  rectangular  region  represented  is  about  5  km  per 
side  with  the  origin  at  the  southwest  corner  of  the  rectangle.  This  region 
can  be  seen  on  the  U.  S.  Geological  Survey  (USGS)  topographic  map, 
Pahramp  quadrangle.  The  origin  is  at  approximately  115  deg  53  min  W, 

36  deg  12  min  N.  Figure  9  shows  the  boundary  of  the  data  base  area. 


Figure  8.  Nonlinear  Texture 


Pahramp  Quadrangle,  Nevada  Series,  USGS 


SECTION  III 


SYSTEM  CONCEPT 


In  this  section,  the  algorithms  selected  for  the  real-time  system  feasibility 
study  will  be  described.  Alternative  algorithms  which  were  considered  are 
also  discussed,  along  with  the  rationale  for  the  selection  which  was  made. 
Functional  requirements  for  the  system  approach  will  also  be  described. 

QUAD  TREE  DATA  STRUCTURE 

The  basic  approach  is  an  outgrowth  of  two  ideas.  One  is  tc  use  a  A -buffer 
to  assemble  components  of  a  scene.  The  second  is  to  use  a  tree  structure 
to  represent  the  terrain  data;  the  root  node  of  the  tree  represents  a  square 
containing  the  entire  geographical  area  for  the  mission.  Each  node  below- 
it  in  the  tree  represents  a  smaller  area  and  thus  a  higher  level  of  detail. 
Figure  10  illustrates  this  quad  tree  structure  and  the  manner  in  which  it 
represents  increasing  level  of  detail  in  the  ground  plane. 

The  algorithm  involved  in  creating  a  picture  is  a  tree  scan.  It  is  not  a  tree 
search  for  some  particular  predetermined  elements  but  a  scan  of  as  many 
branches  as  are  needed  to  find  those  nodes  w-hich  represent  each  pixel  on 
the  display  screen.  Of  course,  a  data  structure  containing  all  possible 
pixels  is  much  too  large  to  implement,  so  the  lowest  nodes  in  the  tree,  the 
terminal  nodes  or  leaves,  are  representations  of  curved  surface  patches. 

The  ground  area  represented  by  each  terminal  node  is  determined  b>  the 
resolution  inherent  in  the  terrain  elevation  data  base.  For  the  scenes 
depicted  in  Figures  7  and  8,  each  terminal  node  represents  a  ground  patch 
of  approximately  100  meters  on  a  side.  Parents  of  the  terminal  nodes 
represent  patches  of  four  times  as  much  area,  the  next  level  of  nodes  16 
times  and  so  on  up  the  tree.  Thus,  each  node  represents  a  quarter  of 
the  area  in  the  ground  plane  represented  by  its  parent  node. 

Each  node  in  the  tree  carries  the  following  data  about  the  surface  area  it 
represents : 

•  Flag  to  indicate  terminal  or  non-terminal  node 

•  X  and  Y  coordinates  of  a  specified  corner  of  the  patch 

•  Length  of  the  sides  of  the  patch 
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Figure  10.  Quad  Tree  Data  Structure 
Non-terminal  nodes  also  carry 

•  Elevations  at  each  of  the  four  corners  of  the  patch 

•  Intensity  (average  of  the  intensities  of  the  four  subnodes) 

•  Pointers  to  its  four  subnodes 
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Terminal  nodes  carry 


•  A  16-entry  array  which  is  a  bicubic  representation  of  the 
surface 

k  Ic 

•  A  texture  map  for  the  surface.  The  map  is  an  2  x  2  array 
of  M  bit  intensities  (for  Figures  7  and  8,  k  =  3  and  M  =  3) 

The  tree  is  scanned  in  depth  first  order.  That  is,  the  algorithm  steps 
through  the  three  from  node  to  node,  beginning  at  the  top,  or  root  node. 
Nodes  near  the  top  represent  large  ground  areas  whose  projections  on  the 
display  screen  will  depend  on  distance  from  the  viewpoint.  Projections 
of  nearby  nodes  will  cover  many  pixels,  so  the  tree  must  be  scanned  to 
greater  depth  for  nearby  nodes.  Distant  nodes  will  project  to  smaller 
picture  areas  and  the  scan  depth  will  be  less.  This  feature  provides  an 
automatic  level  of  detail  selection  as  well  as  minimizing  the  computational 
effort  in  producing  a  picture. 

The  procedure  for  generating  a  picture,  then,  is  as  follows.  Starting  with 
the  top  node,  the  picture  area  corresponding  to  the  node  is  determined  by 
projecting  the  four  corners  of  the  patch  represented  by  the  node  to  display 
screen  coordinates.  A  decision  is  made  to  terminate  the  scan  at  that  node 
or  proceed  to  the  next  level.  Termination  can  occur  under  the  following 
conditions : 

1.  The  projected  area  represented  by  the  node  falls  outside  the 
display  window. 

2.  The  projected  area  faces  away  from  the  viewpoint, 

3.  The  projected  area  is  small  enough  to  display. 

The  result  of  any  of  these  conditions  is  to  terminate  the  scan  at  that  node 
and  discard  all  branches  and  nodes  below  the  termination  node.  When  the 
entire  tree  has  been  scanned  in  this  manner,  the  result  is  the  identification 
of  all  potentially  visible  nodes,  each  node  representing  a  pixel-sized  area 
on  the  display  screen,  and  the  corresponding  display  coordinates.  However, 
some  of  the  potentially  visible  nodes  will  be  occluded  by  other  nodes.  The 
priority,  or  which  nodes  are  actually  visible,  can  easily  be  determined  by 
comparing  the  distance  from  the  viewpoint  to  each  node  where  priority 
conflicts  occur. 


PRIORITY  CONFLICT  RESOLUTION 


The  resolution  of  priority  conflicts  uses  a  combination  of  a  frame  buffer 
and  depth,  or  distance,  buffer  called  a  Z-buffer.  The  frame  buffer  holds 
the  intensities  or  tones  of  display  points  as  they  are  generated.  The  depth 
buffer  stores  the  distances  of  these  same  points  and  resolves  priority 
conflicts  by  keeping  only  the  nearest  point  (both  intensity  and  distance)  when 
conflicts  occur.  An  added  benefit  of  this  approach  is  that  when  the  frame  is 
displayed,  the  distance  to  each  pixel  is  available  for  intensity  modification 
due  to  haze,  fog,  or  other  atmospheric  attenuation. 

SCENE  GENERATION 

The  first  step  required  to  produce  a  picture  or  a  sequence  of  pictures  of 
the  terrain  is  data  base  preparation.  This  step  consists  of  building  the 
data  tree  structure  with  the  data  elements  described  above.  This  step  may 
be  considered  to  be  a  restructuring  of  the  DMA  terrain  elevation  samples, 
augmented  by  texture  data.  The  texture  samples  should  conform  to  the 
DMA  cultural  file  to  represent  the  various  terrain  materials  and  their 
boundaries.  Since  the  scene  generation  algorithm  under  unvestigation  in 
this  study  does  not  depend  on  how  the  texture  arrays  are  actually  generated, 
details  will  not  be  described  here.  The  textures  in  Figures  7  and  8  were 
produced  by  digitizing  aerial  photographs  of  the  terrain  being  represented. 
This  technique  is  only  one  of  several  tin  t  might  be  used. 

The  second  step  is  to  produce  a  skeleton  tree  for  the  present  field  of  view, 
depending  on  viewpoint  location,  direction  and  angular  extent  of  the  field 
of  view.  The  skeleton  tree  is  a  thinned  out  tree  that  contains  pointers  into 
the  full  mission  data  tree  but  has  been  pruned  to  include  only  surfaces  in 
the  current  field  of  view.  This  step  is  repeated  for  each  new  frame.  It 
is  not  necessary,  but  can  substantially  reduce  the  size  of  the  tree  to  be 
scanned  if  the  field  of  view  is  much  smaller  than  the  total  mission  area. 

The  pruning  is  accomplished  simply  by  projecting  the  field  of  view  onto  the 
ground  plane  and  discarding  those  nodes  which  are  completely  outside  the 
field  of  view  boundaries. 

The  terrain  surface  gets  processed  a  node  at  a  time,  starting  with  the  root 
node  of  the  skeleton  tree.  The  depth  buffer  is  initialized  to  a  large  value. 
The  following  steps  get  repeated  until  no  nodes  remain  to  be  examined. 
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1.  Fetch  the  next  node.  If  none  remain,  the  picture  is  ready  for 
display.  Exit. 

2.  Project  the  corners  of  the  node  to  viewpoint  coordinates. 

3.  Count  the  number  of  display  points  covered  by  the  projected  node. 

4.  If  this  is  a  terminal  node  and  there  is  more  than  one  point 
covered,  subdivide  the  node  (Figures  11,  12)  and  go  to  step  1. 

5.  If  there  are  none,  stop  searching  this  branch  of  the  tree  and 
go  to  step  1. 

6.  If  there  is  only  one  and  it  is  behind  the  viewpoint,  discard  the 
node,  stop  searching  this  branch  and  go  to  step  1. 

7.  Compare  the  node  distance  to  that  currently  stored  in  the  depth 
buffer  at  the  covered  display  point.  If  greater,  discard  the  node, 
stop  searching  this  branch  and  go  to  step  1. 


Figure  11.  Ground  Plane  Subdivision  Process 
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8.  Compute  the  displayed  intensity  for  this  node,  store  it  in  the 
frame  buffer,  and  store  the  node  distance  in  the  depth  buffer. 

9.  Go  to  step  1. 

This  process  will  execute  a  depth-first  search  of  the  skeleton  tree  and 
assures  the  examination  of  every  relevant  node.  Any  visible  surface  in 
the  field  of  view  gets  written  to  the  frame  buffer. 

ALGORITHM  DETAILS 

This  subsection  describes  the  algorithms  used  in  Steps  2,  3,  4  and  8. 
Projection  of  Node  Corners 

The  corners  of  a  node  are  known  in  ground  coordinates  as  obtained  from 
the  DMA  terrain  elevation  file.  To  determine  the  position  of  a  patch  on  the 
display  the  four  corner  coordinates  are  translated  and  rotated  into  the 
viewer's  frame  of  reference  and  projected  to  display  coordinates. 

Count  Display  Points  Covered 

To  know  whether  to  stop  searching  some  branch  of  the  tree,  the  number  of 
pixels  covered  by  the  projected  patch  must  be  determined.  If  two  or  more 
pixels  are  covered,  the  node  must  be  subdivided.  The  projected  patch  is 
approximated  with  a  quadrilateral  whose  vertices  are  the  four  projected 
corners  (Figure  13b).  The  number  of  pixel  points  within  this  quadrilateral 
is  then  counted.  Another  method  is  to  count  the  grid  points  in  a  "box" 
determined  by  the  vertices  of  the  quadrilateral  (Figure  13c).  These  methods 
of  determining  the  size  of  the  projected  area  are  evaluated  in  the  tradeoff 
studies  described  in  Section  IV. 

For  non-terminal  nodes,  subdivision  occurs  naturally  as  a  result  of  contin¬ 
uing  the  depth  first  tree  scan.  For  terminal  nodes,  a  bicubic  representation 
of  the  patch  is  maintained  known  as  a  register  square.  Register  square 
subdivision  is  described  in  detail  in  Appendix  R.  Briefly,  the  bicubic 
representation  of  the  four  quarters  of  a  whole  patch  can  be  generated  from 
the  register  square  of  their  parent  with  a  small  number  of  shift  and  add 
operations.  Each  quarter  can  then  form  a  new  register  square  for  further 
subdivision  if  necessary.  Textural  information  is  subdivided  in  parallel 
with  evaluation. 
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Figure  12.  Subdivif 
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Figure  13.  Alternatives  for  Grid  Test 


The  point  intensity  calculation  uses  some  preset  constants,  the  textural 
value  of  the  subdivided  patch,  the  distance  of  the  patch  from  the  viewer, 
and  three  vectors:  the  patch  surface  normal,  the  line  of  sight  from  the 
viewer  to  the  patch,  and,  for  visual  scenes,  a  ray  from  the  sun  to  the 
patch.  Haze  is  an  exponential  function  of  distance.  The  intensities  of  the 
back  sides  of  illuminated  areas  are  thresholded  with  an  ambient  intensity 
constant.  The  point  intensity  computation  is  as  follows: 


T  -  Patch  texture  value 

P 

T  -  Texture  weighting  factor 

w 

1^  -  Ambient  intensity 

Z  -  Distance  of  patch  from  viewpoint 


H  -  Haze  decay  factor 


-  Intensity  of  haze 

S^I  -  Surface  normal  vector 
lAs  -  Dine  of  sight  vector 
SUN  -  Sun  vector 
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R  is  the  dot  product  of  the  surface  normal  vector  and  the  line  of  sight  vector 
times  the  dot  product  of  the  surface  normal  vector  and  the  sun  angle;  this 
product  is  thresholded  by  the  ambient  intensity. 


R  =  max  (I  SN  •  LOS  *  SN  •  SUN) 

T  =  T  *  R  *  T  +  R  *  (1. 0  -  T  ) 
w  p  w 

I  =  point  intensity  =  T  *  e  +  I  *  (1  -  e 

H 

Texture  Generation  and  Representation 

The  DMA  database  specifies  13  classes  of  texture  (surface  materials)  in 
various  bounded  areas  (or  as  a  percentage  of  cover).  The  generated  texture 
information  must  be  added  to  enhance  the  data  base.  Texture  functions  can 
be  generated  either  by  actual  aerial  photographs  of  a  given  texture,  actual 
terrain  data  or  by  computer  synthesis. 

Texture  Function  from  Aerial  Photographs --The  texture  function  to  modulate 
the  diffuse  component  of  the  intensity  model  may  be  generated  with  aerial 
photographs  of  the  texture.  It  may  work  well  for  planar  texture  only  (when 
only  the  "pigment"  of  the  surface  is  textured)  as  when  representing  grass, 
gravel,  etc.  However,  the  digitized  photograph  of  a  three  dimensional 
(nonplanar)  texture  (bumpy  terrain  or  trees,  for  example),  fails  to  simulate 
the  appearance  of  real  texture,  because  the  viewing  angle  and  light  source 
direction  of  the  original  photograph  is  very  rarely  the  same  as  that  of  the 
real-time  simulation. 

Texture  Function  from  Actual  Elevation  Data- -Texture  functions  for  modu¬ 
lating  the  surface  normal  can  be  derived  from  physical  elevation  data  just 
as  the  aerial  photos  generated  the  pigment  function  above.  Examples  are 
terrain  wrinkles,  bumpy  surfaces,  trees  in  clumps,  etc.  This  terrain 
detail  may  or  may  not  have  any  correspondence  with  the  real  terrain  texture 
in  a  given  gaming  area.  The  surface  normal  function  can  be  computed 
directly  from  elevation  data. 


Synthetic  Texture  Function  Generation — Synthetic  texture  functions  have 
the  potential  of  saving  storage  (if  real-time  synthesis  is  achieved)  of  the 
texture  function.  Texture  function  synthesis  can  be  either  deterministic 
or  probabilistic.  A  trivial  example  of  deterministic  synthetic  texture  is  a 
"cultural"  texture  generated  by  an  orderly  (predetermined)  array  of  blocks 
representing  buildings.  This  can  be  either  planar  (painted  on  the  surface) 
or  three  dimensional  and  used  as  the  surface  normal  texture  as  well.  The 
problem  with  this  is  the  need  to  store  a  large  number  of  different  texture 
maps- -otherwise  all  "towns "  would  look  identical. 

Probabilistic  texture  synthesis  allows  the  random  generation  of  texture 
subject  to  a  given  set  of  constraints.  Simplistically,  in  the  above  example, 
the  houses  may  be  required  to  be  in  blocks  but  their  spacing  and  size  may 
be  chosen  from  a  random  number  generator. 

Autoregressive  Texture  Synthesis --Probabilistic  texture  synthesis  works 
for  natural  textures  as  well  (perhaps  better  than  in  cultural  textures).  This 
is  based  on  the  fact  that  texture  functions  can  sometimes  be  modeled  by  an 
autoregressive  function  driven  by  white  noise.  Honeywell  has  investigated 
this  Markovian  texture  synthesis  under  its  ongoing  IR&D  effort.  The  block 
diagram  in  Figure  14  illustrates  this  technique.  The  coefficients  of  the 
autoregressive  function  determine  the  texture  pattern.  The  coefficients 
can  be  determined  from  a  sample  of  the  texture  being  simulated  (say  from 
a  piece  of  an  aerial  photo).  If  the  texture  can  be  adequately  modeled  by  the 
autoregressive  function,  the  synthesized  texture  can  look  remarkably  like 
the  original  texture.  An  additional  bonus- -the  synthesized  texture  can  be 
of  arbitrary  extent  without  repeating  itself  (as  in  certain  "tiling"  schemes). 

Tremendous  data  compaction  can  be  achieved  if  texture  is  represented  and 
synthesized  this  way  in  real  time.  The  only  storage  required  is  for  three 
coefficients  for  a  first  order  two-dimensional  generating  function,  and  the 
initial  values  for  the  filter. 

ALTERNATIVES  AND  TRADEOFFS 

Choice  of  Subdivision  Coordinate  Svstems 

-  —  -  -  -  -  .  » 

The  principal  choice  to  be  made  is  whether  to  subdivide  the  terrain  patches 
in  the  ground-based  coordinate  system  or  in  the  moving  e  -ordinate  system 
centered  at  the  viewpoint,  1.  e.  ,  in  object  space  or  image  space.  If  the 
subdivision  is  done  in  the  ground-based  coordinates,  the  resulting  subpatch 
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a)  Estimation  of  the  Model  Parameters 
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b)  Synthetic  Texture  Generation  Process 

Figure  14.  Markovian  Texture  Synthesis 

corners  must  be  projected  at  each  node  during  the  tree  search.  An  advantage 
to  this  approach  is  that  the  register  squares  can  be  generated  off-line  and 
stored  as  part  of  the  data  base. 

The  image  space  approach  requires  that  the  register  squares  be  computed 
for  every  fra  ne  in  real  time.  Only  the  patch  corners  are  projected, 
resulting  in  far  fewer  multiplication  operations.  For  this  reason,  the  image 
space  subdivision  approach  was  selected.  In  image  space,  the  terrain  is  a 
three-dimensional  vector  parametric  surface,  so  that  three  functions  must 
be  subdivided  simultaneously.  For  the  ground-based  approach,  only  scalar 
subdivision  is  necessary. 

Another  interesting  variation  would  be  to  perform  the  subdivisions  off-line 
in  the  ground-based  coordinate  system.  The  tradeoff  is  clearly  one  of 
increased  memory  requirements  for  off-line  subdivision  versus  increased 
computation  for  on-line  subdivision.  Since  the  ground-based  subdivision 
was  not  selected  for  the  real  time  feasibility  study,  this  tradeoff  was  not 
performed. 
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Binary  Versus  Quaternary  Subdivision 


Quaternary  subdivision  is  the  four-way  subdivision  shown  in  Figure  11. 

An  alternative  is  to  divide  the  patch  into  two  subpatches,  rather  than  four, 
which  can  be  called  "binary"  subdivision.  After  projection  to  the  display 
coordinates,  a  patch  may  be  elongated  since  patches  become  more  fore¬ 
shortened  with  distance  from  the  observer  in  one  direction  (along  the  LOS) 
than  the  other  (perpendicular  to  the  LOS).  In  that  case,  fewer  subdivisions 
are  needed  if  the  patch  is  always  divided  along  its  longest  projected  dimen¬ 
sion.  This  tradeoff  was  done  and  results  are  described  in  Section  IV. 


SECTION  IV 


SYSTEM  DESIGN 


REAL-TIME  SYSTEM  APPROACH 

The  primary  objective  of  the  study  was  to  investigate  the  feasibility  of 
implementing  the  bicubic  spline  subdivision  approach  in  real  time  for 
simulators  and  similar  applications.  This  w'as  done  by  performing  a 
detailed  requirements  analysis  and  system  design  at  the  functional  level. 

Figure  15  shows  the  basic  functions  required  in  the  real-time  system.  Data 
for  the  simulated  mission,  consisting  of  terrain  elevations  and  ground  plane 
texture  are  stored  on  some  mass  storage  device  such  as  disc.  As  the  air¬ 
craft  moves  over  the  terrain,  data  for  the  patches  within  the  field  of  view 
(FOV)  of  the  aircraft  are  transferred  to  the  FOV  memory.  The  FOV 
memory  can  be  visualized  as  a  window  which  moves  across  the  terrain 
with  the  aircraft.  The  size  of  the  window  is  determined  by  the  range  of 
maximum  visibility. 

The  FOV  memory  is  organized  as  a  toroid  with  a  relatively  slow  input  port 
and  a  high  speed  output  port.  For  each  frame,  the  contents  of  the  FOV 
memory  are  scanned  and  processed  to  generate  the  displayed  image. 

The  first  processing  steps  are  referred  to  as  the  Register  Square  Processor. 
This  processor  provides  the  following  functions: 

•  Transforms  terrain  data  from  geocentric  to  eye-centered 
coordinates  with  x  and  y  axes  parallel  to  the  display  x  and  y  axes 

•  Projects  terrain  samples  onto  image  plane 

•  Determines  which  patches  require  subdivision 

•  Computes  bicubic  register  square  values  for  each  patch 
in  the  FOV 

The  register  square  values  for  each  patch  are  the  data  necessary  for  patch 
subdivision.  These  ai  e  stored  in  a  buffer  memory  which  is  organized  as 
a  stack. 
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The  contents  of  the  stack  are  subdivided  by  the  Patch  Subdivider.  Each 
patch  subdivision  results  in  four  subpatches.  Each  subpatch  has  the  same 
data  structure  as  the  parent  patch.  The  subpatches  are  tested  to  see  if 
they  can  be  written  into  the  Z-buffer  by  the  Display  Processor.  Those 
patches  which  must  be  further  subdivided  are  added  to  the  stack. 

The  Display  Processor  determines  the  intensity  or  grey  level  and  distance 
for  each  displayed  point  and  writes  them  into  the  Z-buffer.  The  Z-buffer 
consists  of  a  frame  buffer  which  stores  the  pixel  intensity  values  and  a 
depth  buffer  which  contains  the  distance  to  each  displayed  point.  When 
priority  conflicts  arise,  the  distances  are  compared  by  the  Display 
Processor  and  the  closest  point  is  written  into  the  Z-buffer.  The  frame 
buffer  is  finally  scanned  in  raster  fashion  to  display  the  picture. 

REQUIREMENT  ANALYSIS 

In  order  to  specify  the  computations  and  memory  required  for  each  of  the 
functions  shown  in  Figure  15,  it  was  necessary  to  perform  a  detailed 
requirements  analysis. 

An  overview  of  the  approach  to  this  task  is  shown  in  Figure  16.  The  general 
philosophy  is  to  characterize  the  computational  requirements  of  each  of  the 
algorithm  stages  analytically  as  functions  of  the  flight,  sensor  and  terrain. 
This  characterization  serves  several  purposes: 

•  It  helps  understand  explicitly  how  each  parameter  affects  the 
computational  requirements  of  each  stage  of  the  algorithm. 

•  It  enables  the  prediction  of  the  performance  requirements  for 
varying  flight  sensor  and  terrain  parameters;  the  analysis  is 
not  restricted  to  a  given  sensor  and  mission. 

•  The  computer  simulation  does  not  have  to  be  run  for  all 
possible  cases.  This  is  prohibitively  expensive  besides 
yielding  little  insight  into  the  behavior  of  the  system  under 
varying  conditions. 

The  approach  has  been  to  divide  the  requirements  analysis  problem  into 
two  parts--the  geometry-dependent  part  and  the  algorithm -dependent  part. 
The  former  is  amenable  to  mathematical  analysis  and  allows  direct  use  of 
all  the  flight,  sensor,  and  terrain  parameters.  The  algorithm -dependent 
part  is  harder  to  quantify  analytically,  but  a  few  runs  of  the  simulation  on 
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(Ktpni  >  <>f  Detailed  Requirements  Specifications  Task 


a  limited  test  data  set  have  provided  statistics  on  the  algorithm  performance 
for  each  algorithm  option.  The  results  of  these  statistics  have  been 
combined  with  the  analytically  derived  equations  for  the  geometry-dependent 
part  to  yield  the  detailed  requirements  for  each  algorithm.  These  detailed 
requirements  are  explicitly  parameterized  in  terms  of  the  flight,  sensor, 
and  terrain  parameters. 


Two  alternate  approaches  to  the  perspective  mapping  function  (Figure  17) 
were  considered.  They  are  projection  before  subdivision  and  projection 
after  subdivision.  The  computational  requirements  trade-off  between  the 
two  approaches  suggests  the  former  over  the  latter,  as  discussed  in 
Section  III.  The  detailed  requirements  were,  therefore,  derived  for  the 
projection  before  subdivision  approach. 

An  alternative  approach  to  patch  subdivision--adaptive  binary  subdivision-- 
was  devised,  as  described  in  Section  III.  A  simpler  "box"  test  was  also 
investigated  as  a  grid  coverage  test  in  lieu  of  the  quadrilateral  coverage 
test.  These  have  been  incorporated  into  the  simulation  software.  Three 
versions  of  the  computer  simulations  now  exist:  (a)  binary  subdivision 
with  quadrilateral  test,  (b)  binary  subdivision  with  box  test,  and  (c) 
quaternary  subdivision  with  quadrilateral  test,  the  version  on  STARS. 

The  processing  and  computation  required  to  display  a  patch  is  a  function  of 
the  projected  area  of  the  patch  on  the  screen.  Analytical  equations  were 
derived  (see  Appendix  C)  for  the  projected  area  of  a  patch  as  a  function  of 
the  patch  distance,  the  observer  height  from  the  patch  and  a  measure  of 
the  terrain  "roughness.  "  The  analysis  indicates  that  the  processing  required 
to  display  rough  terrain  can  be  nontrivially  greater  than  that  for  flat  terrain. 
The  equations  form  the  basis  of  the  requirements  analysis  and  were  verified 
by  simulation. 


The  processing  required  to  display  a  patch  as  a  function  of  its  projected 
area  on  the  screen  was  quantified  for  each  of  the  three  versions  of  the 
terrain  representation  simulation.  To  this  end,  the  three  versions  were 
run  with  a  sequence  of  patches  at  varying  distances  and  orientations  from 
the  observer.  Several  statistics  were  gathered  and  analyzed  through 
histograms  and  scatter  plots.  In  particular,  the  orientation  dependence 
of  the  patch  subdivision  process  (binary  vs.  quaternary),  the  relative 
efficiency  of  the  box  test  vs.  the  quadrilateral  test,  the  number  of  sub- 
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Figure  17.  Two  Options  for  Perspective  Rendering 


divisions,  coverage  tests,  intensity  computations,  and  frame  buffer  accesses 
per  patch  as  a  function  of  the  projected  area  of  the  patch  were  quantified. 

A  terrain  roughness  model  was  developed.  This  model  directly  applies  to 
the  non-linear  terrain  representation  approach.  In  this  model,  terrain 
roughness  is  quantified  by  two  parameters:  (a)  the  average  slope  of  a 
terrain  patch,  and  (b)  the  standard  deviation  of  the  elevations  of  the  patch 
corner  points.  The  roughness  parameters  have  been  explicitly  incorporated 
into  the  projected  area  equations.  A  computer  program  was  written  to 
evaluate  the  terrain  roughness  statistics  from  sections  of  the  DMA  data 
base  in  order  to  quantify  the  roughness  parameter  values  for  representative 
terrain  types. 

An  equation  was  also  derived  for  the  distribution  patch  distances.  In 
combination  with  the  equation  for  projected  areas,  the  distance  equation 
provides  the  total  projected  area  as  a  function  of  the  visual  or  sensor  FOV 
and  the  distance  to  the  farthest  patch  to  be  mapped.  These  equations, 
together  with  the  simulation  runs  giving  the  statistics  of  the  algorithm 
performance,  predict  the  number  of  times  an  operation  has  to  be  performed 
at  each  stage  of  the  process  as  a  function  of  sensor,  flight,  and  terrain 
parameters.  The  critical  parameters  are: 

•  Number  of  patches  in  the  FOV 

•  Number  of  patches  which  need  bicubic  subdivision 

•  Number  of  bicubic  subdivisions 

•  Number  of  Z-buffer  accesses 


Number  of  Patches  in  FOV 


Number  of  patches  in  the  FOV  is  given  by  the  equation: 


P  = 


s2> 


2h 


*H{2 

2h2 


where  if  is  the  horizontal  FOV  (in  radians)  of  the  sensor  being  simulated; 
f-2  is  the  distance  to  the  "horizon,  "  the  farthest  patch  being  mapped;  and 
h  is  the  patch  size  on  the  ground. 
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Using  a  nominal  =  60°,  =  15  miles  and  h  =  300  feet,  the  result  is 

P  =  36,  000  patches  in  the  FOV. 

Number  of  Patches  Needing  Bicubic  Subdivision 

The  idea  here  is  that  not  all  patches  which  require  subdivision  have  to  be 
subdivided  using  the  bicubic  subdivision.  If  the  projected  patch  area  is 
small  (for  example,  <m,  a  threshold),  then  the  register  square  computation 
can  be  eliminated,  and  processing  proceeds  directly  to  the  patch  subdivision. 
The  subdivision  on  these  patches  would  then  be  simply  linear  interpolation. 
Thus,  the  expensive  bicubic  fits  and  register  square  computation  on  these 
patches  are  not  required.  The  number  of  patches  requiring  bicubic  sub¬ 
division  is  given  by: 

B  =  —  9  m 

m  mTT  t 

H 

Here  m  is  the  area  of  a  patch  in  display  grid  points  beyond  which  it  requires 
bicubic  subdivision;  0m  is  the  terrain  roughness  parameter  (the  average 
patch  slope),  and  n2  is  the  total  number  of  raster  points. 

The  patch  tilt  0  and  the  patch  altitude  were  histogrammed  for  two  sections 
of  the  DMA  data  base  (Nevada).  One  was  a  large  100  km  x  100  km  section 
and  the  other  was  the  mountainous  5  km  x  5  km  section  that  was  used  for 
image  generation  examples.  These  histograms  are  shown  in  Figures  18a, 

18b,  19a,  and  19b.  It  can  be  seen  that  the  average  0,  0m,  for  the  mountainous 
terrain  is  about  10°;  this  is  the  value  used  in  the  worst  case  example  for  the 
requirements  analysis. 

Assuming  0m  =  10°  (rough  terrain),  n  =  1000  (high  resolution  display),  and 
m  =  10,  Bm  11,000.  However,  for  the  hardware  specification  it  was 
assumed  that  all  patches  in  the  60°  FOV  (~40,  000)  require  bicubic  sub¬ 
division. 

Total  Number  of  Patch  Subdivision,  Ns 


The  number  of  patch  subdivisions  required  is  proportional  to  the  total 
projected  area.  A,  of  the  patches.  Each  of  the  three  versions  of  the  basic 
algorithm  has  a  different  proportionality  constant  which  was  determined  by 
simulation.  The  three  versions  are: 
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stogram  of  Patch  Altitudes 
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•  Binary  subdivision  with  quadrilateral  grid  test 

•  Binary  subdivision  with  box  test 

•  Quaternary  subdivision  with  quadrilateral  grid  test 

Each  of  the  three  versions  was  simulated  with  a  set  of  patches  at  different 
distances.  The  number  of  patch  subdivisions  and  Z -buffer  accesses  and  the 
projected  patch  area  were  gathered  for  each  patch.  These  statistics  were 
plotted  in  scatter  plots  and  analyzed  for  each  algorithm  version  and  patch 
orientation  with  respect  to  the  viewpoint  (Figures  20,  21,  and  22).  From 
these  the  constants  of  proportionality  that  link  the  total  number  of  subdivi¬ 
sions  and  Z-buffer  accesses  for  each  algorithm  to  the  total  projected  area 
(which  is  an  explicit  function  of  the  sensor  flight  and  terrain  parameters) 
were  determined.  It  is  observed  that  the  box  test  for  grid  coverage  appears 
to  require  many  more  subdivisions  than  other  versions  (using  quadrilateral 
test  for  the  45°  patch  orientation  case). 

Assuming  that  all  patch  orientations  are  equally  likely,  the  results  were: 

{1.  5 A  (Quaternary) 

1. 3A  (Binary  with  quadrilateral  test) 

2.  2  A  (Binary  with  box  test) 

Note  that  these  are  equivalent  two-way  subdivisions  (one  quaternary  sub¬ 
division  is  computationally  equivalent  to  tw'o  binary  subdivisions).  The  total 
projected  area  is  derived  in  Appendix  C  and  is  given  by: 
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where  '  is  the  height  of  the  viewpoint.  Under  simplifying  approximations, 
this  becomes  (see  Eq.  (2-15,  Appendix  C): 

.  2  f  12  6m] 

A  --  n  1  +  —  — 


raster  points,  which  explicitly  shows  the  effect  of  6 
parameter. 
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,  the  terrain  roughness 


PROJECTED  AREA  (PIXELS) 
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NUMBER  OF  SUBDIVISIONS/UNIT  AREA 


O  PATCH  ANGLE  =  45° 
O  PATCH  ANGLE  =  90° 


Figure  20.  Scatter  Plot  of  Projected  Area--Binary 
(Quadrilateral  Test) 
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PROJECTED  AREA  (PIXELS) 
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Figure  21.  Scatter  Plot  of  Projected  Area--Binary  (Box  Test) 
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PROJECTED  AREA  (PIXELS) 
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O  PATCH  ANGLE  =  45° 

O  PATCH  ANGLE  =  90° 

Figure  22.  Scatter  Plot  of  Projected  Area --Quaternary 
(Quadrilateral  Test) 


Ih 


r 


Assuming  9m  =  15°  for  a  worst  case  terrain  roughness,  a  vertical  FOV 
jr  =  30°  gives  a  projected  area: 

A  =  n2  [2.27] 

which  implies  that  a  rough  terrain,  on  the  average,  can  have  more  than  two 
hidden  surfaces  which  map  to  the  same  point  on  the  screen.  This  result 
and  conclusion  are  discussed  in  detail  in  Appendix  C. 

As  a  baseline,  the  quaternary  subdivision  with  the  quadrilateral  test  was 
chosen  for  further  hardware  specification.  Although  binary  subdivision 
required  fewer  subdivisions,  the  resulting  data  structure  would  be  more 
complicated.  The  number  of  quaternary  subdivisions  for  the  given  para¬ 
meters  is  then  given  by: 

N  =  1. 5  x  (1/2)*  x  2.  27  n2 
s 


or 

g 

N  a  1.7x10  /frame 
s 

g 

for  a  display  resolution  of  10  raster  points. 

Figures  23a,  23b,  23c,  and  23d  show  the  predicted  projected  area  vs. 
distance  curves  for  various  heights  and  average  patch  tilts  0m.  Figures 
24a,  24b,  and  24c  show  the  ratio  of  the  average  predicted  area  to  the 
average  measured  area  from  the  simulation  for  each  patch  distance,  height, 
and  tilt  8m.  Note  that  the  predicted  results  track  the  measured  results 
closely.  This  is  important  because  the  requirements  analysis  is  founded  on 
the  projected  area  equations. 

Number  of  Z -Buffer  Accesses 


The  number  of  Z-buffer  accesses  is  also  a  function  of  the  total  projected 
area.  In  fact,  for  both  versions  of  the  patch  subdivision  algorithm  with 
the  exception  of  the  box  test,  this  constant  is  unity  (that  is,  one  Z-buffer 
access  /'unit  projected  area).  This  has  been  verified  by  computer  simulation. 
Hence, 


The  factor  1/2  arises  because  this  is  the  number  of  four-way  subdivisions. 
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Figure  23a.  Projected  Area  vs.  I)istance--o 
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PROJECTED  PATCH  AREA  IN  PIXELS 
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Figure  24a.  Comparison  of  Predicted  and  Measured 
Projected  ^reas--0  = 
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Figure  24c.  Comparison  of  Predicted  and  Measured 
Projected  Areas--S  =  15J 


NZ  = 


A  for  all  algorithm  versions,  except, 

1. 72 A  for  box  test  (0*  patch  orientation) 


For  the  rough  terrain  parameter  assumed, 

NZ  =*•  2.27  x  106. 

REAL-TIME  IMPLEMENTATION  REQUIREMENTS 


Figure  25  shows  the  storage  and  throughput  requirements  for  the  various 
functional  units  of  the  system,  based  on  the  analysis  described  in  the 
preceding  section.  It  should  be  noted  that  these  requirements  are  conser¬ 
vative,  based  on  worst  case  requirements  in  most  instances. 

Implementation  requirements  were  specified  for  each  functional  block  of 
Figure  25.  The  approach  taken  was  to  assume  fully  parallel  computation 
and  to  use  pipeline  architectures  where  possible  for  the  individual  blocks 
and  from  block-to-block.  In  a  fully  pipelined  system,  the  throughput  is 
determined  by  the  slowest  part  of  the  pipeline. 


Field-of-View  (FOV)  Memor\ 


The  FOV  Memory  provides  temporary  storage  for  terrain  elevation  samples 
for  all  potentially  visible  patches.  All  patches  for  the  mission  are  stored 
in  a  mass  memory  such  as  a  magnetic  disc.  However,  access  to  mass 
memories  is  too  slow  for  the  pipelined  Register  Square  Processor,  so  the 
FOV  Memory  is  necessary.  The  contents  of  the  FOV  Memory  are  continually 
updated  as  the  aircraft  flies  across  the  terrain.  A  maximum  of  about  500 
patches  must  be  added  and  subtracted  for  each  new  frame.  Approximately 
one  terrain  elevation  sample  of  16  bits  is  required  for  each  patch. 


The  FOV  memory  is  organized  like  a  toroid  in  that  addresses  which  exceed 
the  physical  size  of  the  memory  wrap  around  and  begin  again  at  the  initial 
addresses.  This  organization  allows  for  rapid  scanning  of  all  patches  within 
the  FOV  and  easy  updating  with  new  patches  as  the  aircraft  moves  across 
the  terrain. 


Once  every  frame,  the  terrain  elevation  data  are  scanned  out  of  the  FOV 
Memory  and  processed  by  the  Register  Square  Processor. 
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Register  Square  Processor 


The  Register  Square  Processor  subsystem  is  a  pipelined  processor  which 
converts  all  relevant  terrain  patch  data  into  bicubic  form  in  image  screen 
coordinates  for  subdivision  and  display. 

Hie  separate  functional  units  of  the  Register  Square  Pipeline  Processor  are 
shown  in  Figure  26.  The  first  block  (Stage  1)  transforms  the  terrain  ele¬ 
vation  samples  to  perspective  space.  The  required  calculations  are  shown 
in  Figure  27.  A  pipelined  hardware  implementation  of  this  function  is 
shown  in  Figure  28. 

The  pipeline  contains  nine  multipliers,  nine  adders,  four  registers,  and  two 
dividers  and  requires  five  clock  cycles  for  execution.  The  throughput  is 
determined  by  the  slowest  operation,  the  divide.  Assuming  a  conservative 
16-bit  divide  time  of  200  nsec  (typical  current  generation  off-the-shelf 
device  speeds),  all  40,000  terrain  samples  in  a  60°  FOV  can  be  projected 
in  8  x  10-3  sec.  Inputs  to  the  projection  hardware,  in  addition  to  the  terrain 
samples,  are  aircraft  position  coordinates,  heading,  roll,  and  pitch.  The 
outputs  are  the  terrain  elevation  samples  expressed  in  image  screen 
coordinates. 

To  determine  the  bicubic  patch  coefficients,  it  is  necessary  to  collect  the 
terrain  elevation  samples  in  overlapping  groups  of  16  samples /patch.  This 
is  the  function  of  the  next  stage,  called  the  circular  buffer  (Stage  2).  For 
each  patch,  the  16  sample  vectors  defining  the  patch  are  used  in  Stage  3  to 
determine  the  bicubic  fit  coefficients,  which  in  turn  are  used  in  Stage  4  to 
compute  the  "register  squares.  "  The  register  squares  are  the  quantities 
used  in  the  patch  subdivision. 


FOO* 

F«V  "EMORY 


STAr.E  1  S’AGE  2  STAGE  1  S’AGf  i 


Figure  26.  Functional  Elements  of  Register  Square  Pipeline  Processor 


vare  Block  Diagram 
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The  computations  required  for  the  bicubic  fit  coefficients  are  described  in 
Appendix  A  and  are  summarized  in  Figure  29.  This  stage  requires  the 
multiplication  of  the  4x4  array  of  16  projected  terrain  sample  vectors 
surrounding  each  patch  by  a  constant  4x4  matrix  M  (premultiplication) 
and  its  transpose  (post  multiplication). 

A  pipelined  hardware  approach  to  the  bicubic  fit  (Stage  3)  is  shown  in  Figure 
30.  Each  x,  y,  and  z  component  requires  32  identical  modules,  each 
consisting  of  four  multipliers  and  three  adders.  Sixteen  registers  are  also 
required  for  the  4x4  constant  matrix  M.  Six  clock  periods  are  required 
for  completion  of  each  patch.  Again,  the  throughput  rate  is  conservatively 
estimated  at  5  x  10®  patches /second. 

The  output  of  the  bicubic  fit  stage  consists  of  a  4  x  4  matrix  A.  There  is  an 
A  matrix  for  each  of  the  three  components.  The  elements  of  the  A  matrix 
are  denoted  a. .. 

The  next  step  (Stage  4)  consists  of  computing  the  16  register  square  entries 
from  the  a^j  for  each  of  the  three  components.  The  definitions  of  the 
register  square  entries  are  shown  in  Figure  31  and  are  derived  in  Appendix 
B.  A  pipelined  implementation  of  this  stage  is  shown  in  Figure  32.  A  total 
of  27  adders  and  49  registers  are  required  for  parallel  implementation. 

Each  patch  requires  four  clock  periods. 

In  summary,  the  Register  Square  Processor  performs  the  functions  of 
perspective  projection  and  the  fitting  of  the  bicubic  polynomial  to  the  terrain 
samples.  The  output  is  in  the  form  of  three  arrays  of  16  numbers  which  are 
used  directly  by  the  succeeding  bicubic  subdivision  algorithm.  This  pro¬ 
cessor  is  fully  pipelined  and  can  process  an  FOV  of  40,000  patches  in 
8  msec.  The  hardware  required  consists  of  393  multipliers,  two  dividers, 

327  adders,  and  72  registers. 

Register  Square  Stack 

The  outputs  from  the  Register  Square  Processor,  as  well  as  the  previously 
subdivided  patches  which  require  further  subdivision,  are  stored  in  a  memory 
called  the  Register  Square  Stack  until  they  can  be  processed  by  the  Patch 
Subdivider.  This  memory  can  be  organized  as  either  a  first  in-first  out 
(FIFO)  or  last  in-first  out  (LIFO)  stack. 
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Figure  29,  Bicubic  Fit 


Figure  30.  Pipelined  Bicubic  Fit  Block  Diagram 
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Figure  31.  Register  Squares  Prior  to  Subdivision 

It  is  difficult  at  this  time  to  predict  the  size  of  the  stack  needed.  It  will  be 
a  function  of  the  input  rate  of  the  patches,  but,  more  importantly,  it  will  be 
the  rate  at  which  subdivided  patches  are  returned  for  further  subdivision. 

The  total  memory  capacity  required  is  given  by 

.  ,  16  16  3  number  of  patches  in  the 

(Register  X  (Bits)  X  (x,y,  z)  X  stack  at  a  given  time 
Squares) 

In  the  worst  case,  assume  that  one  patch  fills  the  FOV.  Tbe  worst  case 
enumerates  the  subdivision  tree  breadth  first.  Since  the  10®  terminal 
patches  do  not  go  back  to  the  stack,  all  nodes  of  the  quaternary  tree  have 
to  be  stored.  Thus,  4/3  x  10® /4  =  0.  33  x  10®  subpatches  may  be  in  the 
stack.  Hence,  the  maximum  size  is  0.33  x  10®  x  768  =  256  megabits. 

Using  64Kbit  random  access  memory  (RAM)  chips  devices,  this  memory 
requires  4000  integrated  circuits.  These  can  be  arranged  into  blocks  to 
reduce  the  input-output  (I/O)  speed  requirement  of  each  stage.  Note  that 
this  is  the  upper  limit.  It  is  quite  possible  that  the  memory  requirement 
can  be  substantially  less  than  this.  Tliis  can  be  verified  by  Monte  Carlo 
simulations. 
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PATC  DIVIDER 


The  algor  n  'r  tch  subdivision  is  described  in  Appendix  B  and  is  shown 

in  Fign  ■■  .  ^o.t  hat  subdivision  produces  nine  2x2  register  squares 

from  tm  .^ur  2  x  °  iste.  squares  of  the  parent  patch.  The  nine  new 
register  squares  -  uped  to  provide  four  new  sets  of  register  squares 

correspc  ding  to  the  r  .  subpatches. 


Figures  34a  and  34b  show  the  parallel  hardware  implementation  of  the  sub¬ 
divide  algorithm.  It  requires  93  adders  and  318  registers.  With  this  degree 
of  parallelism,  four  clock  periods  are  required  for  subdivision  of  a  bicubic 
patch  x(u,  v),  y(u,  v),  z(u,  v).  For  a  clock  period  of  100  nsecs,  three  patch 
subdivision  pipelines  are  required. 


Tester 


After  subdivision,  each  subpatch  is  tested  to  see  whether  it  should  be  further 
subdivided  or  passed  on  to  the  Display  Processor.  The  present  test  criterion 
is  that  the  patch  in  question  should  surround  one,  and  only  one,  display  grid- 
point.  Patches  which  surround  more  than  one  gridpoint  are  returned  for 
further  subdivision. 

Two  alternatives  to  approximating  the  projected  patch  for  the  grid  test  have 
been  considered:  the  quadrilateral  and  box  approximations  (Figure  13). 

The  box  approximation,  as  can  be  seen,  tends  to  overestimate  the  number 
of  grid  points  surrounded  by  the  patch  but  results  in  a  simpler  implemen¬ 
tation.  Figures  35  and  36  show  the  implementation  of  the  box  grid  test. 

For  maximum  throughput,  four  grid  testers  are  required  for  each  Patch 
Subdivider,  for  a  total  of  12. 

Display  Processor  and  Z -Buffer 

The  Display  Processor  is  a  pipelined  processor  that  determines  the  tone 
or  intensity  for  each  display  subpatch  and  assembles  the  picture.  In  the 
picture  assembly  process,  any  priority  conflicts  are  resolved  point-by- 
point  using  the  Z-buffer. 
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NOTE: 


SUBSCRIPT  I  PINT  If  Its  Out 
OF  FOUR  NCW  REGISTER 
SQUARE  VALUES  IN  ONE 
Of  NINE  NEW  REGISTER 
SQUARES,  DENOTED  BV  CIRCLED 
NUMBER. 


Figure  34a.  Subdivide  Pipeline  Part  I 
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Figure  34b.  Subdivide  Pipeline  Part  II 
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Figure  35.  Grid  Point  Test 


Figure  36.  Max/Min  Circuit  as  Part  of  Grid  Point  Test 


The  Z-buffer  is  a  memory  containing  the  distance  to  each  displayed  subpatch 
as  well  as  its  tone  or  intensity.  Priority  conflicts  are  resolved  by  comparing 
the  distances  to  the  two  subpatches  in  question.  The  subpatch  which  is  closest 
is  entered  into  the  Z-buffer  and  the  other  is  discarded.  Two  Z -buffers 
allow  readout  of  one  frame  to  the  display  while  the  next  frame  is  being 
assembled. 

Thirty-two  parallel  channels  are  required  for  the  Z-buffer.  This  estimate 
is  based  on  the  estimated  requirement  for  approximately  2.  3  x  10®  intensity 
computations  every  frame.  This  results  in  an  intensity  computation  every 
14.5  nsec. 

Assume  that  dense  MOS  memories  will  be  used  for  these  memory-intensive 
stages,  with  180  nsec  access  times  and  300  nsec  cycle  times.  Consequently, 
the  Z-buffer  and  intensity  buffer  updates  will  take 

1 80  nsec  Read 
40  nsec  Compare 
300  nsec  Write 
520  nsec  Total 

The  update  speed  required  is  14.  5  nsec /update.  Thus,  the  number  of 
memory  channels  needed  is 

520  nsec 

-  =-  36 

14.5  nsec 

Each  of  the  36  channels  will  address  a  specific  sector  of  the  display  screen. 
Pixels  will  ue  assigned  to  one  of  the  36  channe's  on  the  basis  of  x,  y,  screen 
coordinates. 

The  intensity  computation  is  a  function  of  the  angles  between  the  LOS,  the 
sun  and  the  patch  normal,  and  the  texture  function  and  an  atmospheric 
attenuation  term  based  on  patch  distance.  The  texture  function  is  deter¬ 
mined  by  table  look-up,  using  the  patch  ground  coordinates  for  addressing. 

The  intensity  computation  pipeline  is  shown  in  Figure  37.  This  processor 
requires  29  multipliers,  18  adders,  76  registers,  and  two  table  look-up 
read-onlv  memories  (PROMS)  per  channel. 
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Figure 


Only  seven  channels  of  intensity  computations  are  required.  The  intensity 
computation  is  pipelined  and  the  pipeline  clock  rate  is  10  MHz,  so  that  the 
effective  throughput  rate/channel  is  10^. 

The  required  throughput  pixel  rate  is  about  2.  3  x  3  x  10  /second,  resulting 
in  a  need  for  seven  parallel  channels.  This  is  the  product  of  the  number  of 
pixels/frame  (10®),  the  number  of  frames /second  (30)  and  the  additional 
factor  (2.3)  to  account  for  extra  computations  from  conflicts  (Appendix  C). 

The  total  parts  count  for  the  system,  excluding  the  host  computer  and  mass 
memory,  is  contained  in  Table  1.  This  count  also  does  not  include  control 
circuitry,  multiplexers,  etc.,  whose  specification  requires  a  more  detailed 
design.  Allowing  1000  devices  for  these  additional  functions,  the  total  parts 
count  is  approximately  9000.  Of  this  total,  approximately  two-thirds 
consist  of  memory  devices,  including  4000  64K  x  1  devices  for  the  Register 
Square  Stack.  As  discussed  previously,  this  estimate  is  probably 
conservative,  but  improvement  requires  further  analysis  by  Monte  Carlo 
simulation. 

TIMING  ANALYSIS 

Figure  25  showed  the  throughput  requirements  in  each  stage  of  the  pipelined 
processor  in  terms  of  the  number  of  principal  operations  for  each  stage. 
These  requirements  are  reflected  in  the  subsequent  estimation  of  the  through¬ 
put  and  hence  hardware  requirements  for  each  pipelined  stage.  The  impli¬ 
cation  is  that  the  preliminary  architecture  discussed  supports  the  system 
throughputs  shown  in  Figure  25.  This  subsection  addresses  the  total 
transport  delay  in  the  pipelined  cascaded  stages. 

The  total  transport  delay  is  reckoned  from  the  real  time  FOV  memory  to 
display.  This  measures  the  time  elapsed  from  a  new  platform  position 
input  (from  the  flight  simulator)  to  the  display  of  this  new  position  on  the 
screen. 

The  transport  delay  T^  is  the  sum  of  the  individual  pipelined  stage  delays: 

T  =  T  +T  +T  +T  +T 

D  REG  STACK  SUB  TEST  DISP 
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PARTS  COUNT  FOR  PIPELINED  APPROACH 


where 


and 


T 

REG 

tstack 

T 

SUB 

T 

TEST 

T 

DISP 


delay  through  register  square  computation, 

delay  through  the  stack, 

delay  through  the  patch  subdivision, 

delay  through  the  tester, 

delay  through  the  display  processor 


Each  of  these  delays  is  estimated  below. 

T 

REG:  Referring  to  Figure  26,  this  includes  projection,  circular  buffer, 
bicubic  fit  and  register  square  computation.  The  projection  (Figure  28) 
requires  five  basic  stages  (subtract,  multiply,  two  adds  and  a  divide)  of 
delay.  Assuming  200  nsec  of  delay  per  stage  (corresponding  to  the  worst 
stage--the  divide),  this  implies  200  x  5  =  1  usee.  The  circular  buffer  adds 
a  further  delay  of  four  rows  of  patches:  200  x  4  x  200  nsec  =  160  usees. 
Here  we  have  assumed  200  nsec /patch  and  200  patches /row  (corresponding 
to  40,000  patches  in  the  FOV).  The  bicubic  fit  has  a  delay  of  six  stages 
(Figure  30),  or  6  x  200  nsec  =1.2  usee.  The  register  square  pipeline 
(Figure  32)  has  a  delay  of  four  stages  (800  nsec). 

Hence  T_„  =  1  +  1 60  +  1 .  2  +  0.  8  =  163  u,secs. 

R 

T 

ST AC.K:  This  is  basically  the  delay  associated  with  the  buffering  of  Ihe 
patches,  and  is  variable,  depending  on  the  size  of  the  patches  received 
(range  to  the  patch).  The  worst  case  delay  occurs  when  the  stack  fills 
up  with  patches  and  the  succeeding  pipelines  are  busy  and  cannot  accept 
patches  from  the  stack.  This  upperbound  is  1  frame  time  (33  msec),  i.  e.  , 

T  =  33  msec  (max) 

AL  K 

T 

■SUB :  The  subdivide  pipeline  entails  a  total  delay  of  eight  stages  (see 
Figures  34a  and  34b),  or  1.6  ..secs  at  200  nsec/stage 

TS!  B  =  1  •  6  ^secs 
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T 

TEST:  This  consists  of  the  gridpoint  test  (Figures  35  and  36)  or  eight 
stages  at  200  nsecs /stage 

Hence  TTEST  =  1.6  usees 

T 

DISP:  The  intensity  computation  pipeline  involves  (Figure  37)  17  stages 
at  200  nsecs.  or 

tdisp  =  3'4>lSecs 

Hence,  the  total  worst  case  transport  delay 

=  163  +  33,000  +  1.6  +  1.  6  +  3.4 
D 

=  33,  169.  6  nsecs 

Rounding  off  to  the  nearest  larger  frame  time,  the  worst  case  transport 
delay 


T 

=  2  frame  times 

D 
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SECTION  V 


CONCLUSIONS  AND  RECOMMENDATIONS 


The  principal  conclusion  of  this  study  is  that  real-time  implementation  of 
the  approach  is  feasible  for  aircraft  visual  and  sensor  simulations.  The 
hardware  requirements  are  not  unreasonable,  although  it  seems  that  they 
could  be  significantly  reduced  through  additional  iterations  of  the  design 
tradeoffs.  In  particular,  the  Register  Square  Stack  member  which  was 
conservatively  sized  at  256  megabits,  could  be  substantially  reduced. 

Improvements  could  be  made  in  the  test  to  determine  whether  a  subpatch 
should  be  displayed.  The  number  of  patch  subdivisions  could  be  reduced  to 
the  point  where  perhaps  only  one  Patch  Subdivider  pipeline  would  suffice  if 
subdivision  was  not  required  to  single  picture  element  size. 

Since  the  number  of  patches  to  be  processed  for  fixed  patch  size  increases 
as  the  square  of  the  distance,  it  would  be  desirable  to  use  variable  patch 
sizes.  Distant  patches  need  not  be  as  small  as  nearby  patches  where  greater 
detail  is  necessary.  This  would  not  only  provide  variable  resolution  as  a 
function  of  viewing  distance  but  would  reduce  the  size  of  the  FOV  Memory. 
Together  with  the  previous  suggestions  for  improving  the  grid  test,  this 
would  also  reduce  the  size  of  the  Register  Square  Stack.  The  primary 
obstacle  is  the  possibility  of  edge  effects  and  artifacts  where  patches  of 
different  sizes  are  juxtaposed  in  the  image  plane. 

Some  hardware  reduction  may  be  possible  in  the  Register  Square  Processor. 
With  a  fully  pipelined  design,  the  duty  cycle  of  this  stage  is  only  25  percent. 
Since  it  basically  consists  of  a  sequence  of  4  x  4  matrix  multiplications, 
some  time-sharing  of  hardware  may  be  possible.  This  stage  is  also  a  prime 
candidate  for  exploiting  very  large  scale  integrated  circuit  (VLSIC) 
technology. 

Further  work  should  be  done  on  the  Display  Processor.  The  computation 
of  tone  or  intensity,  including  the  representation  of  texture,  is  by  no  means 
optimized  at  this  point. 
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Finally,  more  work  is  required  on  anti-aliasing  with  this  approach.  Some 
anti-aliasing  is  inherent  in  the  averaging  of  texture,  based  on  projected 
patch  size.  However,  this  approach  still  exhibits  severe  aliasing  problems 
along  the  horizon  and  along  any  boundary  which  separates  a  cultural  object 
from  a  terrain  surface  patch. 

An  anti-aliasing  approach  has  been  under  investigation  at  Honeywell  for 
anti-aliasing  with  the  bicubic  subdivision  algorithm.  First,  note  that  area 
sampling,  i.e. ,  averaging  of  patch  contributions,  is  needed  only  at  an 
object's  silhouette  boundaries  where  multiple  surface  contributions  have 
to  be  averaged.  This  is  because  texture  aliasing  can  be  eliminated  by 
prefiltering  the  texture  pattern  on  a  patch  depending  on  its  distance  from 
the  observer  (farther  patches  are  filtered  more  than  the  nearer  ones). 
Therefore,  point  sampling  the  filtered  texture  can  be  used  everywhere 
except  at  object  or  terrain  silhouettes.  In  fact,  this  texture  filtering  is 
already  performed  implicitly  in  the  current  implementation  of  the  Honeywell 
terrain  representation  program  on  STARS  at  AFHRL  (by  having  several 
levels  of  texture  detail  and  choosing  the  level  corresponding  to  the  patch 
distance).  Hence,  texture  aliasing  can  be  eliminated  by  adaptive  texture 
filtering  in  real  time  before  the  patches  are  subdivided  and  tested. 

Further,  subpatches  which  come  from  silhouettes  can  be  recognized  by 
their  surface  normals  (which  are  orthogonal  to  the  line  of  sight).  Hence, 
pixels  which  correspond  to  the  silhouette  edges  can  be  treated  separately 
from  all  the  rest,  and  multiple  surface  averaging  applied  only  at  these 
pixels.  Because  the  silhouette  pixels  comprise  only  a  small  fraction  of 
the  image,  this  insight  results  in  a  greatly  simplified  anti-aliasing  algorithm 
for  object  space  (depth  buffer  algorithms). 

As  explained  above,  texture  prefiltering  on  each  patch  eliminates  the  need 
for  area  sampling  everywhere  but  at  the  silhouette  pixels.  Hence,  the 
unmodified  point  sampling  version  of  the  algorithm  is  applied  to  all  patches 
(and  subpatches)  which  are  not  tangential.  However,  tangential  (silhouette) 
subpatches  are  treated  differently.  If  a  tangential  subpatch  is  mapped  to  a 
certain  pixel,  the  Z-buffer  entry  corresponding  to  the  pixel  is  not  compared 
or  replaced  with  the  new  subpatch.  Instead,  a  tangential  subpatch  arriving 
at  a  pixel  signals  that  the  point  is  a  potential  silhouette  point.  The  current 
content  of  the  Z-buffer  is  replaced  with  a  pointer  to  a  linked  list  which  will 
contain  information  about  the  current  and  future  contributions  to  that  pixel. 

The  current  contents  of  the  intensity  buffer  and  Z-buffer  for  this  pixel  and 
the  information  from  the  new  (silhouette)  patch  form  as  the  first  two  elements 
of  this  list.  At  the  end  of  the  first  pass,  all  nonsilhouette  points  are  completely 
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determined.  Moreover,  the  information  from  all  surfaces--the  area 
(A),  depth  (Z)  and  intensity  (I)  is  available  in  the  linked  lists  for  each 
potential  silhouette  pixel.  This  is  then  used  to  perform  averaging  and 
priority  determination  for  the  silhouette  pixels. 

The  above  approach  satisfies  the  following  goals  for  efficiency: 

1.  The  expensive  multi-surface  area  sampling  is  done  only  for 
the  potential  silhouette  points,  eliminating  the  jagged  edges. 

2.  The  use  of  linked  lists  for  storing  information  on  the  sil¬ 
houette  points  eliminates  the  need  for  multiple  frame  and 
Z -buffers,  and  uses  only  as  much  memory  as  needed  to 
store  the  silhouette  information. 

3.  Texture  prefiltering  enables  the  computationally  simple 
(conventional)  point  sampling  algorithm  to  be  used  almost 
everywhere  in  the  image,  without  texture  aliasing. 
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APPENDIX  A 


BICUBIC  SPLINE  INTERPOLATION 


Currently,  polygons  are  used  to  represent  terrain  elevation  in  computer- 
image  generation.  Approximation  of  terrain  with  polygons  produces  a 
faceted  effect  and  a  silhouette  made  of  straight- line-segments.  Curved 
surface  segments  or  "patches"  can  be  used  instead  of  polygons.  If  the 
patches  can  be  joined  together  with  slope  continuity,  a  picture  of  the 
surface  can  be  made  continuous  both  in  shading  and  silhouette. 


The  most  widely  used  nonlinear  representation  is  the  bicubic  spline,  an 
extension  of  the  cubic  spline  interpolation  for  functions  of  one  variable. 

In  general,  a  spline  curve  is  a  piecewise  analytical  function  whose  pieces 
are  chained  together  with  continuity  requirements  on  the  first  few  deriva¬ 
tives  imposed  at  the  joints.  The  pieces  are  basic  curve  shapes  with 
adjustable  parameters  for  each  piece.  These  parameters  are  chosen  to 
satisfy  the  imposed  continuity  requirements.  Polynomial  curves  are 
usually  used  for  the  pieces.  The  cubic  polynomial  f(u)  =  auJ  +  bu2  +  cu  +  d 
is  the  lowest  degree  polynomial  with  a  sufficient  number  of  parameters  to 
provide  continuous  differentiability. 


A  curve  in  space  can  be  represented  in  parametric  form  by  the  vector 
equation 


r(u) 


x(u) 

y(u) 

z(u) 


where  each  of  the  components  can  be  represented  by  a  cubic  defined  for 
the  parameter  u  in  the  interval.  Usually,  for  convenience,  0  £  u  s  1.  In 
matrix  notation,  each  of  the  cubics  can  be  expressed  as 
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A  surface  patch  is  a  function  of  two  variables  u  and  v. 


F  (u.  v) 


x(u,  v) 
v(u,  v) 
z(u,  v) 


Terrain  elevation  data  is  single- valued  and  is  provided  in  the  DMA  data 
base  as  z(x,y)  over  a  (nearly)  uniform  mesh  X:,  Vj.  By  associating  x  -*  u 
and  y  -»  v,  the  parametric  representation  of  onl\  the  vertical  component 
z(u,  v)  is  ne  eded. 


The  matrix  notation  for  z(u,  v)  is 


z(u,  v) 


n  o 

u  u"  u  1 


11 

ai2 

a  1  3 

"14 

21 

a22 

a2  3 

a24 

31 

a32 

"l3  ' 

a34 

41 

%2 

a43 

a44 

-  u  A  v 


where  the  a-  -  are  coefficients  of  the  equation  just  as  a.  b,  c,  and  cl  w  ere 
in  the  univariate  case.  The  problem  is  to  find  the  coefficients. 

Following  Catmull^,  only  local  methods  will  be  considered,  i.e.  ,  each 
datum  only  affects  the  coefficients  of  nearbv  patches. 


One  method,  using  an  array  of  16  samples  of  z(u,  v).  might  simultaneously 
solve  the  resulting  linear  system  of  equations  in  the  unknown  a- This 
method  would  result  in  a  bicubic  patch  which  would  pass  through  all  the 
data  points  and  interpolate  the  surface-  h  ctween.  In  addition  to  being 
computationally  undesirable,  this  method  would  result  in  :  surface  that 
would  be  discontinuous  in  the  first  derivative  at  *hc  patch  boundari*  s. 
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Consider  again  the  univariate  equation.  Fiur  data  samples,  say  fj,  f9, 
fg,  f 4,  are  required  to  determine  the  four  coefficients  a,  b,  c,  d.  The 
coefficients  can  be  related  to  the  data  samples  by  some  4x4  matrix  M 


Therefore, 


This  concept  can  be  trivially  extended  to  the  bivariate  case: 


T 

Note  that  A  '»!  Z  M  .  For  the  cons  i  d<  red  earlier  of  fitting  •«  bitub 

through  16  po.nts,  the  matrix  M  is  given  b\ 


This  technique  circumvents  the  need  to  solve  a  system  of  16  linear  equations, 
but  the  problem  of  discontinuity  in  the  derivative  of  the  surface  at  the  patch 
boundaries  remains. 


THE  B-SPLINE 

The  cubic  B-spline  provides  an  alternative  which  provides  continuity  of  the 
second  derivative.  In  general,  except  when  the  points  are  collinear,  the 
B-spline  does  not  interpolate  the  data  samples  but  rather  approximates 
them.  This  property  provides  some  smoothing  of  the  input  data  samples. 
Consider  the  four  points  f  ,  f  ,  f  ,  f  • 

1  u  J  *1 


A  cubic  curve  segment  can  be  generated  w'hich  does  not  pass  through  any  of 
the  points.  Similarly  another  segment  of  curve  can  be  generated  using  the 
points  fg,  fg,  f^,  and  a  fifth  point  fg  which  connects  to  the  first  segment 
with  second  derivative  cor.tinuity. 


86 


For  the  B-spline  cubic,  the  matrix  M  is  derived  from  a  set  of  basic  functions. 
The  weights  of  the  basic  functions  are  chosen  to  satisfy  the  continuity 
conditions,  resulting  in: 


M 


2 


1_ 

6 


3 

-6 

0 

4 


-3 

3 

3 

1 


1 

0 

0 

0 


There  are  some  interesting  geometric  properties  of  these  piecewise  curves. 
They  lie  everywhere  within  the  convex  hull  of  the  polygon  of  the  defining 
vertices.  At  u  =  0  or  u  =  1,  the  curve  passes  through  a  point  P  which  is  the 
1/3  point  of  the  median  of  the  triangle  formed  by  three  sequential  vertices. 
The  first  derivative  vector  at  P  is: 


P’ 


2 


and  the  second  derivative  vector  at  P  is  P"  =  (f.  -fj+j)  +  (fi+2  " 


i 

j 


i 

i 


i 


f 


When  the  three  vertices  f.,  f^,  ^1+2  are  co^inear»  then  the  triangle  is 
degenerate  with  P  on  the  tine.  The  tangent  vector  P'  and  second  deriva¬ 
tive  P"  are  also  along  the  line. 

B-splines  have  the  "variation  diminishing"  property  which  means  that  they 
approximate  linear  functions  exactly,  as  demonstrated  above,  and  the 
approximation  is  always  "smoother"  than  the  function  it  approximates. 

THE  CATMULL-ROM  CUBIC  SPLINE 


This  spline  interpolates  the  data  points  and  has  continuity  of  first  derivative 
at  the  joints.  Consider  the  four  points  f^,  f^,  f^,  f^: 


0  1  u  -*■ 


A  cubic  can  be  generated  that  interpolates  from  point  f  to  f  Consider  a 
fifth  point  f^. 


Another  segment  of  curve  is  generated  using  f2«  fg.  f<j»  The  two  segments 
join  with  first  derivative  continuity.  For  the  Catmull-Rom  spline,  the 
matrix  M  is 
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APPENDIX  B 

BICUBIC  SPLINE  SUBDIVISION 
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APPENDIX  B 


BICUBIC  SPLINE  SUBDIVISION 


The  binary  bicubic  spline  subdivision  algorithm  described  here  is  an 
adaptation  of  an  algorithm  developed  by  E.  Catmull  in  his  thesis  published 
in  December  1974  entitled,  "A  Subdivision  Algorithm  for  Computer  Display 
of  Curved  Surfaces.  "  Catmull  develops  an  algorithm  for  subdividing  a 
cubic  patch  and  its  first  derivatives  in  two  orthogonal  directions  and 
discusses  computational  requirements  for  the  algorithm. 

CUBIC  SUBDIVISION 

Consider  the  cubic  function 

3  2 

f(T)  =  aT  +bT+cT  +  d 

which  is  valid  for  t-hST^t+h.  By  adding  f(t  +  h),  the  value  of  f  at  the 
center  point  is  given  by 

f(t)  =  [f(t  -  h)  +  f(t  +  h)]/2  =  h2  (3at  +  b)  (B_1} 

=  [f(t  -  h)  +f(t  +  h)]/2  -  g(t) 

The  correction  term 

g(t)  =  h2  (3at  +b)  (B-2) 

can  also  be  expressed  in  terms  of  its  values  at  the  ends  of  the  interval 

g(t)  =  [g(t  -  h)  +  g(t  +  h)]/2  (B-3) 

Also  note  that,  if  h  is  the  width  of  the  n^  interval,  at  the  (n+l)St  binary 

subdivision  h  , ,  =  n  / 2  and 
n+1  n 

gn+l(T)  =  gn(T)/4  (B'4) 
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Equations  (B-l)  through  (B-4)  provide  the  basis  for  the  cubic  spline  sub¬ 
division  algorithm  described  below.  Using  matrix  algebra  notations, 

define 


c 

— n- 


1 


(t) 


f(-r) 

g^jH  (2x1) 


Then  to  scale  the  correction  term  for  the  subdivision  in  progress,  form 


To  compute  the  function  and  correction  term  for  the  center  point  of  the 
interval,  t  s  t  s  t 

i  « 


c 

— n 


B^ause  the  elements  of  H  and  Q  are  powers  of  1/2,  a  binary  search  of  a 

M  IV 

cubic  spline  is  computationally  reduced  to  shifts  and  adds.  This  desirable 


property  carries  over  into  subdivision  of  a  bicubic  spline  function. 
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BICUBIC  SUBDIVISION 

As  described  in  Appendix  A,  the  bicubic  spline  for  a  terrain  patch  can  be 
described  by  the  surface 


where 


z(u,  v)  =  u'  M  Z  M1 


v 

0  s  u  <;  1 


0  <s  v  ^  1 


(B-6) 


(u,  v)  are  independent  variables  (corresponding  to  the  orthogonal  coordinates 
(x,  y)  for  the  bicubic  surface  within  the  patch;  the  prime  denotes  a  trans¬ 
pose.  M  is  a  (4  x  4)  matrix  of  second  derivatives,  and  Z  is  a  (4  x  4)  matrix 
of  surface  elevations  surrounding  and  including  the  patc£T  in  question,  i.  e. , 


Z 

where,  if  1 


=  {z.O  =  {z(x.,  yj} 

£  i  and  i  <.  4,  the  terrain  patch  in  question  is  bounded  by 
<:  x  s  x 

O 

*  y  5  y3- 


To  develop  the  extension  of  the  subdivision  algorithm  to  two  dimensions, 
define 


A  =  M  Z  M1  =  {a..}  (4x4) 

~  U 

and  consider  the  vector 

s'(u)  -  u'  A  s  [SJ-S2  s3  s4] 


with 


4  (4-i) 

s .  =  l  ,  u  a. . 
J  i  =  l  U 


Then 


s(u,  v)  =  ^V  ^Sj 


(B-7) 


93 


For  constant  v,  each  Sj  is  a  cubic  in  u  which  has  a  correction  term  gj 
associated  with  it.  In  particular,  s(u,  vQ)  has  the  correction  term 

4  (4-i) 

g(u.v)  =  rvQ  ij(u)  (b-8) 

j  =  l 

For  constant  u,  both  s(uQ,  v)  and  g(u0>  v)  are  cubics  in  v.  Hence  each  has 
its  own  correction  term,  say  cg(u  ,  v)  for  s  and  cg(uQ,  v)  for  g.  Following 
Catmull,  these  four  functions  can  Be  arranged  in  a  "register  square:" 

v 


Moving  in  the  v  direction,  cs  corrects  s,  and  Cg  corrects  g.  Conversely, 
in  the  u  direction,  g  corrects  s  and  c  corrects  cg.  Based  on  the  cubic 
subdivision  algorithm  described  previously,  these  functions  are  explicitly 


In  particular,  the  register  squares  corresponding  to  the  corner  of  the 
terrain  patch  prior  to  subdivision  (h  =  h^  =  1)  are  given  by: 
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4 

4 

4  4 

4 

r  a 

J-i  4j 

E  a2j 

3*1 

E  T  a.. 

i=j  j=l  1J 

I  (3a,  .  +  a.,.) 

j=l  11  2) 

4 

3<3all  +a2l’ 

+  (3a12  +a22> 

3a41  +  a42 

3a21  +a22 

E  (3an  +ai2) 

1  =  1 

4 

a44 

a24 

E 

i=l  ai4 

3a!4+a24 

4 

a42 

a22 

I 

i=l  ai2 

3a!2  +a22 

With  these  preliminary  definitions,  the  bicubic  subdivision  algorithm  can 
be  expressed  in  matrix  notation  as  follows.  For  any  patch  or  subpatch,  let 
the  register  squares  for  the  corners  of  the  patch  be  described  by  a  set  of 
(2x2)  matrices  in  the  form 


(n  —  1 ) 

r  (n-1) 

01 

~  11 

(n-1) 

r  (n-1) 

00 

~  10 

In  preparation  for  subdivision,  the  correction  terms  in  each  register  square 
are  first  scaled  by  the  operator 


1  0 


0  I 


to  produce  scaled  register  squares 


Subdivision  will  produce  a  (3  x  3)  block  of  (2  x  2)  register  square  matrices 
of  the  form 


(n) 

(n) 

(n) 

£oi 

1 

£n 

(n) 

(n) 

(n) 

c  l 

~o,  * 

c  1  , 

~  2,  2 

£i,t 

(n) 

(n) 

(n) 

in 

o 

o 

£*.0 

£  io 

In  terms  of  the  subdivision  operator. 


Q  =< 


*  -*  4  -* 

o*o  * 


the  new  register  squares  resulting  from  the  bicubic  subdivision  are  then 
given  by 


(n)  (n)  (n) 

~*.o  =  *•  £oo  ~10 


(n) 

cn  i  =  Q 

~0,  f  ~ 


C 

~01 


(n) 

31 

(n) 


o00  ^ 


(n) 


(n)  (n) 


5  S' ;  s*.i  =  fa,!  £n  1  Q’ 


(n) 

£1.*  =2 


^  <np 

£11 

(n) 

r 

~  10 


^  (nP 


<n) 

£*,*  =  f  So.* 


(n)  (n)  . 

£i.*  }  8  S’  =  S 


~*.i 


(n) 

~*,0 
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Any  of  the  four  resulting  distinct  sets  of  four  register  squares  can  then  be 
subdivided  by  repretion  of  the  algorithm. 


It  is  worth  noting  that  the  off-diagonal  correction  terms  in  C  decrease  as 
4~n,  and  the  diagonal  correction  term  decreases  as  16~n.  '’t’his  fact  can 
perhaps  be  used  to  reduce  the  computational  load  by  easing  into  a  linear 
interpolation  (i.  e. ,  stop  updates  of  the  correction  terms)  after  a  few  sub¬ 
divisions  in  the  binary  search  have  been  performed. 

Surface  Derivatives 

In  addition  to  searching  the  surface  of  the  bicubic  patch  for  an  intersection 
with  the  LOS,  the  first  derivatives  of  the  surface  are  also  required  to 
compute  the  direction  cosines  of  the  surface  normal  at  the  point  of  inter¬ 
section  and  thereby  the  reflected  intensity  for  a  given  sun  angle. 

If  one  considers  the  cubic  and  its  derivative 

3  2 

f(-r)  =  a  t  +  b  t  +  c  t  +  d 
2 

f  (t)  =  3a  t  +  2b  t  +  c 

T 

the  derivative  at  the  center  of  the  interval  t-h  <;  t  <  t  +  h  can  be 
expressed  as  either 
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r 


fT(t)  =  ff(t-fh)  =  f(t-h)]  /  (2h)  -  ah2 
or 

f  (t)  =  [f  (t-h)  +  f  (t4h)]  /  2  -  3ah2 

T  T  T 


(B-9) 


Catmull  suggests  the  former  approach,  but  the  latter  is  more  attractive 
for  the  two-dimensional  subdivision  algorithm  because  the  same  subdivision 
algorithm  can  be  used  for  the  surface  derivatives  as  for  the  surface.  More 
importantly,  the  algorithm  based  on  the  first  equation  requires  that  the  four 
patches  bounding  the  patch  in  question  be  subdivided  simultaneously.  For 
example,  along  an  outer  boundary  such  as  u  =  0,  computation  of  su  (0,  v) 
requires  knowledge  of  s(-hn,  v),  which  is  not  contained  in  the  patch  being 
subdivided. 

Proceeding  as  in  the  case  of  the  surface  subdivision,  and  using  Equation 
(B-9)  as  the  basic  one- dimensional  algorithm,  consider  the  first  derivative 
in  the  u  direction: 


s  (u,  v)  =  [3u  2u  1  0]  A  v  = 
u  ~  — 


4 

V 

l-i 

3=1 


(4-j) 

V  J  £ 

i=l 


(4-i)  u 


(3-i) 


a. . 


For  constant  v,  s  (u,  v  )  is  a  c  uadratic  in  u  with  correction  term 
u  o 


g  (u.  v  )  =  3h  T. 
u  o  u  .  , 
3=1 


On  the  other  hand,  for  constant  u,  both  su(uQ,  v)  and  gu(u0,  v)  are  cubics 
in  v  and  have  the  appropriate  cubic  form  for  the  correction  coefficients. 

Based  on  these  observations,  the  initial  set  of  register  squares  for  the  u 
derivative  of  the  terrain  patch  can  be  quickly  determined. 


98 


Register  Squares  for  s  (u,  v)  Patch 
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A  similar  consideration  of  the  v  derivative  of  the  surface,  sv(u,  v),  yields  a 
set  of  register  squares  which  are  similar  to  those  for  su(u,  v)  but  with  the 
(i,  j)  indices  interchanged,  the  C.  .  transposed,  and  the  diagonal  C  inter¬ 
changed,  ~ 


(0.  1) 
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APPENDIX  C 


PROJECTED  AREA  OF  A  PATCH  AT  THE  SCREEN 


The  first  step  is  to  derive  a  key  relationship  between  the  orientation  of  a 
patch,  relative  position  of  the  observer,  and  the  projected  area  of  a  patch 
at  the  screen.  Much  of  the  further  analysis  is  based  on  this  relationship. 

The  imaging  geometry  is  shown  in  Figure  C-l.  The  patch  is  approximately 
at  a  distance  Tlfrom  the  observer  (who  is  in  the  £-n  plane).  The  observer 
is  at  height  £  relative  to  the  patch.  The  normal  to  the  patch  makes  an 
angle  8  with  the  vertical,  and  its  projection  on  the  £Tl  plane  makes  an  angle 
0  with  the  §  axis  (cylindrical  notation).  The  line-of-sight  makes  an  angle  0 
with  the  horizontal.  The  patch  size  is  h  x  h. 

For  a  flat  terrain,  the  projected  area  of  the  patch  in  angular  subtense  at 
the  screen  is  given  by: 


,  h  A  2 

Area  =  cos  A  rad 


(C-l) 


where  A  is  the  angle  made  by  the  line-of-sight  vector  X  with  the  patch 
normal  n. 

Now,  cos  A  =  7  •  n  where: 

X  -  [cos  0,  0,  sin  g] 

-*  T 

and  n  =  [sin  0  cos  0,  sin  6  sin  0,  cos  ©] 


Then: 


cos  A  =  [sin  0  cos  0  +  cos  8  sin  9  cos  0]  . 


(C-2) 


Substituting  (C-2  into  C-l): 


Area  =  [sin  8  cos  9  +  cos  0  sin  0  cos  0] 


(C-3) 
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Projected  Area  of 


For  small  p,  the  observation  angle  can  be  approximated  by: 
sin  p  =-  £/? 

COS  p  =“  1 

Similarly,  sin  0^6  and  cos  9  =*  1  for  small  0. 

Hence  the  area  of  the  patch  becomes: 


Equation  (C-4)  gives  the  projected  area  of  a  patch  at  a  specific  tilt  orienta¬ 
tion  (0,  0),  distance  Tl,  and  relative  height  £  from  the  observer.  In  practice, 
of  course,  patches  can  have  different  tilts  9,  orientations  0,  and  heights  £ 
with  respect  to  the  observer.  Hence,  to  find  the  average  projected  area  of 
a  patch  at  a  given  distance,  statistical  models  of  the  parameters  9,  0, 
and  £  have  to  be  generated.  This  statistical  model  gives  rise  to  a  means 
of  characterizing  the  terrain  roughness. 

TERRAIN  ROUGHNESS  MODEL 


The  terrain  is  composed  of  square  patches  (of  side  h)  each  with  an  angle  9 
to  the  vertical,  an  angle  0  to  the  vertical  plane  containing  the  patch  and  the 
observer,  and  height'  below  the  observer  (Figure  C-l).  Note  that  0  denotes 
whether  a  patch  is  tilted,  and  0  determines  whether  it  is  tilted  toward  or 
away  from  the  observer.  The  following  assumptions  can  be  made  about  the 
distribution  of  9,  0  and  £  over  the  entire  field  of  view: 

•  9,  0  and  £  are  independent  random  variables 

•  9,  0  and  £  are  stationary  processes  in  the  field  of  view--that  is, 
the  distribution  of  the  parameters  is  not  a  function  of  the  location. 


It  can  be  further  assumed  that  0  is  uniformly  distributed  between  0  and  2tt. 
This  means  it  is  equally  likely  that  the  observer  is  looking  at  a  patch  from 
any  direction  of  the  compass.  The  distribution  of  0,  the  tilt  of  a  patch  to 
the  vertical,  need  not  be  completely  specified.  It  is  completely  characterized 
for  purposes  of  this  analysis  by  its  mean  value  9m.  Finally,  £  is  assumed 
to  be  a  normal  random  variable  with  mean  m<  and  standard  deviations 
m,  and  ar  are  measured  over  the  significant  region  of  the  field  of  view. 
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The  terrain  roughness  is  then  completely  characterized  by  9m,  the  average 
slupe  of  a  patch,  and  ct  ,  the  standard  deviation  of  the  patch  elevations  over 
a  small  area.  - 

THE  AVERAGE  AREA  OF  A  PATCH  AT  DISTANCE  X 


From  the  terrain  roughness  model,  the  average  area  of  a  patch  in 
equation  (C-4)  is  computed  by  integrating  over  the  distributions  of  9,  0 
and  Q.  Since  these  parameters  are  assumed  to  be  independent,  the  order 
of  integration  is  immaterial. 


The  result  in  equation  (C-4)  is  first  averaged  over  all  possible  orientations 
of  the  observer  with  respect  to  the  patch.  Since  each  orientation  0  is 
equally  likely,  the  probability  density  of  0  is: 


p(0) 


tt-  0  <  0  <  2tt 
4it 

0  elsewhere 


(C-5) 


Since  the  "back"  of  a  patch  must  be  subdivided  if  it  is  visible,  negative 
areas  (of  a  patch  facing  away  from  observer)  should  be  counted  as  positive 
areas.  Therefore,  the  average  over  0  should  be  written: 


E  farea] 
0 


9  cos  0 


(C-6) 


where  !  •  !  denotes  absolute  value,  and  E  [  ]  denotes  expected  value  or 
average. 


Note  that  as  0  varies  from  0  to  tt,  cos  0  takes  on  values  from  between  +1  and 
-1.  Therefore,  the  above  quantity  is  difficult  to  evaluate  in  closed  form 
unless  the  specific  value  of  £/ti  is  known.  Instead,  Schwartz's  inequality 
gives:  r 


E  [area]  s  ~ 

T) 


r 

A 


n 


+  E 


9  cos  0 


(C-7) 


which  gives  an  upper  bound  (conservative  estimate)  on  the  area. 
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Now: 


cos  <t  d  0 


(C-8) 


Since  this  equation  is  linear  in  0,  averaging  over  all  possible  0  and 
assuming  0  has  a  mean  value  of  6^  gives: 


E 

0 


[area] 


-  e 

n  m 


A 


(C-10) 


Note  that  the  first  term  in  the  above  equation  corresponds  to  the  flat  earth 
case.  The  second  term  modifies  the  result  for  a  rough  terrain. 

DISTRIBUTION  OF  PATCH  PROJECTED  AREAS 


Equation  (C-10)  gave  the  average  projected  area  of  a  patch  at  a  distance  71 
and  relative  height  '  from  the  sensor.  To  compute  the  total  projected  area 
at  the  screen  :!ue  to  all  the  patches  in  the  FOV,  the  distribution  of  the  patch 
distances  must  be  determined--that  is,  how  many  patches  are  at  a  given 
distance  from  the  observer.  The  number  of  patches  in  the  strip  d "H  at  a 
distance  ~|  is  given  by: 

p(-)d^  =  -9-  dT  (C-ll) 

h"  H 

The  total  projected  area  is  obtained  by  integrating  (C-10)  and  (C-ll)  over 
the  near  and  far  distances  in  the  FOV,  that  is: 


1  05 


A<VV 


A(H)p(tl)d  n 


/ 


£  .  t  +  2  e 

n  tt  m  71 
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H 


T] 


l2  "  *1 


v2 


dll 


+  w  6m  *H  Ln  (  T[  '  lrad2 


(C-12) 


Because  ecluation  (C-12)  can  be  written  as: 


A(VV  ■  ^  *H  +  !  6m  *H  L"  (if 


(C-l  3) 


The  projected  area  is,  of  course,  a  function  of  the  near  distance  (the 
distance  to  the  bottom  of  the  FOV).  The  worst  case  is  when  the  horizon  is 
at  the  top  of  the  screen  as  shown  in  Figcre  C-2,  Then  the  vertical  field  of 
view  determines  -tj  because  from  Figure  C-2,  ^  =  £/-£.  . 


Figure  C-2.  Imaging  Geometry  for  the  Worst  Case 
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Then  equation  (C-13)  becomes: 


A 


*v*H 


(C-l  4) 


Assuming  an  nxn  element  raster  over  the  FOV  (f  •  i|r  ),  the  area  in 
number  of  raster  elements  is  given  by: 


TT 


(C-l  5) 


Equation  (C-l 5)  explicitly  shows  the  total  projected  area  for  a  rough  terrain 
when  the  entire  field  of  view  is  covered  with  the  ground  patches.  In  the  flat 
earth  case,  it  is  obvious  that  A  =  n“  which  is  intuitively  correct.  With  rough 
terrain  (0m  >  0),  because  of  hidden  surfaces  more  than  one  point  on  the 
ground  maps  to  the  same  point  on  the  screen.  This  is  the  reason  why  the 
total  projected  area  can  be  greater  than  n2  raster  elements.  To  see  how 
rough  terrain  affects  the  quantity  in  equation  (C-15),  assume  =“  500  m, 

-tg  =  20  km,  and  a  vertical  field  of  view  =  30°.  Table  C-l  shows  the 
total  projected  area  as  a  function  of  6  ,  the  roughness  parameter. 

Table  C-l  shows  that  for  0m  =  15  (extremely  rough  mountain  terrain)  the 
total  projected  area  is  2.  27  times  the  corresponding  result  for  the  flat 
earth  case  because  of  hidden  surfaces. 


TABLE  C-l.  AREA  VS.  8 

m 


e ' 

m 

0 

A/n" 

0 

1 

2 

1.  16 

5 

1.39 

7.  5 

1. 59 

10 

1.78 

12 

1.94 

15 

2.27 
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APPENDIX  D 

NUMBER  OF  PATCHES  NEKDiNO  BICUBIC  SUBDIVISION 
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APPENDIX  D 


NUMBER  OF  PATCHES  NEEDING  BICUBIC  SUBDIVISION 


As  discussed  earlier,  the  number  of  potentially  visible  patches  in  the  field 
of  view  is  very  large  (200,000  to  400,000  to  a  range  of  50  miles).  However, 
not  all  these  patches  have  to  be  subdivided  to  render  them  on  the  display. 
Whether  or  not  a  patch  has  to  be  subdivided  depends  on  its  projected  area 
in  terms  of  the  number  of  screen  grid  crossings  it  covers.  If  it  covers  more 
than  one  screen  grid  intersection,  it  has  to  be  subdivided.  Furthermore,  as 
mentioned  in  the  discussion  of  the  project  and  subdivide  approach,  it  may 
be  desirable  to  use  the  bicubic  subdivision  algorithm  on  only  the  nearer-in 
patches  which  span  a  large  number  of  screen  grid  intersections.  More 
distant  patches  may  be  subdivided  using  a  linear  subdivision.  This  saves 
the  expensive  real-time  bicubic  fit  and  register  square  computation  on  the 
large  number  of  patches  which  can  be  subdivided  using  linear  subdivision. 
Here,  the  relationship  is  quantified  between  the  number  of  patches  which 
require  bicubic  subdivision  and  the  sensor,  terrain,  and  flight  parameters 
which  affect  this  number. 

Assume  that  a  patch  will  be  subdivided  using  bicubic  subdivision  if  its 
projected  area  (in  terms  of  grid  intersections)  is  greater  than  m.  Let  l 
be  the  distance  at  which  a  patch  subtends  on  the  average,  m  picture  grid 
points.  The  number  of  patches  in  the  field  of  'view  up  to  -t  is  given  by  the 
equation  to  be: 


B 

m 


II 


i  2 
111 


( 1 )-  1 ) 


The  screen  grid  intersections  subtended  bv  a  patch  is  ;■  function  of  the 
resolution  of  the  sensor  and  the  tilt  of  the  patch  with  respect  to  the  observer 
and  .the  distance  to  the  observer.  In  Appendix  C,  the  average  projected  area 
(in  angular  subtense)  of  a  patch  was; 


A 


9 

h“ 


TT 


(C-10) 
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Where:  h  is  the  size  of  the  patch; 

T|  is  the  distance  to  the  observer; 

£  is  the  height  of  the  observer  relative  to  the  patch;  and 
0m  is  the  "average"  slope  of  the  patch  with  respect  to  the 
horizontal,  a  function  of  the  terrain  roughness. 

At  large  distances  and  rough  terrain,  the  second  term  in  Equation  (C-10) 
dominates  the  first  term,  that  is : 


The  number  of  screen  grid  intersections  is  given  by: 


(D-2) 


(D-3) 


Where  di|r  is  the  resolution  of  the  simulated  sensor.  The  distance  at  which 
a  patch  subtends  m  screen  grid  intersections  is  given  by  equations  (D-2)  and 
(D-3)  as: 


(D-4) 


From  which: 


(D-5) 


The  number  of  patches  B  in  the  field  of  view  at  a  distance  less  than 
is  given  by  (D-l)  and  (D-4)  to  be: 


B 

m 


e 

m 

nm 


(D-6) 
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If  =  n'  (the  number  of  resolution  elements  in  a  line  of  the  display). 

Equation  (D-6)  can  be  rewritten  as: 


B 

m 


1 

mn 


2 


n 


(D-7) 


The  result  in  Equation  (D-7)  is  interesting  because  it  explicitly  shows  the 
dependence  of  Bm  on  the  surface  roughness  (9m),  horizontal  field  of  view 
(f^j),  and  number  of  display  elements  (n^).  For  a  sensor  with  =  60°, 
and  assuming  a  typical  rough  terrain  of  8  =  10°,  Bm  =“  212000/m  is 

derived  from  Equation  (D-7),  If  the  threshold  m  =  5,  only  21200  patches 
have  to  be  subdivided  using  the  bicubic  subdivision  algorithm;  that  is,  the 
bicubic  fits  and  register  squares  must  be  computed  in  real  time  on  these 
patches.  The  remaining  patches  in  the  field  of  view  can  be  subdivided 
linearly. 


At  this  time,  no  subjective  simulation  results  have  been  obtained  as  to  the 
largest  value  of  m  which  will  not  produce  artifacts  due  to  the  linear  inter¬ 
polation.  A  factor  m  =  5  appears  reasonable,  especially  in  view'  of  the 
conservative  1000  x  1000  display  resolution  assumed. 
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